MatrixTmplDec.hxx

Go to the documentation of this file.

/*
 * Copyright (c) 2001-2004 MUSIC TECHNOLOGY GROUP (MTG)
 *                         UNIVERSITAT POMPEU FABRA
 *
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */


#ifndef _MatrixTmplDec_
#define _MatrixTmplDec_

#include <iosfwd>
#include "Array.hxx"
#include "CLAM_Math.hxx"

namespace CLAM
{

        template <class T>
        class MatrixTmpl
        {
        public:
                MatrixTmpl();
                MatrixTmpl(unsigned int dim1, unsigned int dim2);
                MatrixTmpl(const MatrixTmpl<T>& originalMatrix);
                ~MatrixTmpl();
                int GetNumRows() const {return mNumRows;}
                int GetNumColumns() const {return mNumColumns;}
                int GetNumElements() const {return mNumRows*mNumColumns;}
                const Array <T>& GetBuffer() const {return MatrixBuffer();}
                Array <T>& GetBuffer() {return MatrixBuffer();}

                T Sum() const
                {
                        T sum=0;
                        for (unsigned int i=0; i<(mNumRows * mNumColumns); i++)
                                sum += MatrixBuffer()[i];
                        return sum;
                }

                T Max() const
                {
                        T max = MatrixBuffer()[0];
                        for (unsigned int i=1; i<(mNumRows*mNumColumns); i++) 
                                if (MatrixBuffer()[i] > max)
                                        max = MatrixBuffer()[i];
                        return max;
                }

                T Min() const
                {
                        T min = MatrixBuffer()[0];
                        for (unsigned int i=1; i<(mNumRows*mNumColumns); i++) 
                                if (MatrixBuffer()[i] < min)
                                        min = MatrixBuffer()[i];
                        return min;
                }

                float Mean() const {return (float)(Sum())/(mNumRows*mNumColumns);}

                // Print the matrix
                void Print() const;

         // Determinant
                float GetDet() const
                {
                        CLAM_ASSERT(mNumRows == mNumColumns,
                                "Determinant only valid for square matrices");
                        float sum = 0;

                        if (mNumColumns == 1)
                                return (*this)(0,0);

                        for (unsigned int i=0; i< mNumColumns;i++) // For the first row
                        {
                                sum += MatrixBuffer()[i] * pow((TData)-1,(TData)i) * (GetSubmatrix(0,i).GetDet());
                        }
                        // i+1 : the matrix index are from [0..N-1] 
                        return sum;
                }       

                // Get Transposed Matrix
                MatrixTmpl<T> GetTrans()
                {
                        MatrixTmpl<T> ret(mNumColumns, mNumRows);
                        for (unsigned int i=0; i<mNumColumns; i++)
                                for(unsigned int j=0; j<mNumRows; j++)
                                        ret(i,j) = (*this)(j,i);
                        return ret;
                }

                // Transpose the matrix
                void Trans()
                {
                        MatrixTmpl<T> ret(mNumColumns, mNumRows);
                        for (unsigned int i=0; i<mNumColumns; i++)
                                for(unsigned int j=0; j<mNumRows; j++)
                                        ret(i,j) = (*this)(j,i);
                        (*this) = ret;
                }
                // Reset the matrix (set all the elements to 0)
                void Reset()
                {
                        MatrixTmpl<T> ret(mNumRows, mNumColumns);
                        for (unsigned int i=0; i<mNumRows; i++)
                                for(unsigned int j=0; j<mNumColumns; j++)
                        ret(i,j) = 0;
                        (*this) = ret;
                }


                // Invert the matrix
                void Invert()
                {
                        CLAM_ASSERT(mNumRows == mNumColumns,
                                "Inverse can only be calculated for square matrix");
                        MatrixTmpl<T> ret(mNumRows, mNumColumns);
                        MatrixTmpl<T> aux(mNumRows-1, mNumColumns-1);
                        for (unsigned int i=0; i<mNumRows; i++)
                                for (unsigned int j=0; j<mNumRows; j++)
                                {
                                        aux = GetSubmatrix(i,j);
                                        ret(i,j) = pow((TData)-1,(TData)i+j) * aux.GetDet();
                                }
                        ret = (1 / GetDet()) * (ret.GetTrans());
                        (*this) = ret;
                }

                // Get the inverse of a matrix
                MatrixTmpl<T> GetInverse() const
                {
                        CLAM_ASSERT(mNumRows == mNumColumns,
                                "Inverse can only be calculated for square matrix");
                        MatrixTmpl<T> ret(mNumRows, mNumColumns);
                        for (unsigned i=0; i<mNumRows; i++)
                                for (unsigned j=0; j<mNumColumns; j++)
                                        ret(i,j) = pow((TData)-1,(TData)i+j) * ( GetSubmatrix(i,j).GetDet() );
                        ret = (1/GetDet()) * (ret.GetTrans());
                        return ret; 
                }

                // Delete a Row and return a new Matrix Oject with reduced dims
                MatrixTmpl<T> GetDelRow(unsigned int row) const
                {
                        MatrixTmpl<T> ret (mNumRows-1,mNumColumns);

                        for (unsigned i=0;i<row;i++) // copy unshifted part
                        {
                                for(unsigned j=0;j<mNumColumns;j++)
                                        ret(i,j) = MatrixBuffer()[i*mNumColumns+j];
                        }
                        for (unsigned i=row;i<mNumRows-1;i++) // shift rest one row up
                        {
                                for(unsigned j=0;j<mNumColumns;j++)
                                        ret(i,j) = MatrixBuffer()[(i+1)*mNumColumns+j];
                        }
                        return ret; 
                }

                // Delete Column
                MatrixTmpl<T> GetDelColumn(unsigned int column) const
                {
                        MatrixTmpl<T> ret (mNumRows,mNumColumns-1);

                        for (unsigned i=0;i<mNumRows;i++) // shift matrix one column left
                        {
                                for(unsigned j=0;j<column;j++) 
                                        ret(i,j) = MatrixBuffer()[i*mNumColumns+j];
                        }

                        for (unsigned i=0;i<mNumRows;i++) // shift matrix one column left
                        {
                                for(unsigned j=column;j<mNumColumns-1;j++)
                                        ret(i,j) = MatrixBuffer()[i*mNumColumns+j+1];
                        }
                        return ret;
                }

                // Get a Submatrix(i,j) deleting the i row and the j column
                MatrixTmpl<T> GetSubmatrix(unsigned int i, unsigned int j) const
                {
                        MatrixTmpl<T> aux( GetDelRow(i) );
                        MatrixTmpl<T> ret( aux.GetDelColumn(j) );
                        return ret;
                }

                // Convert matrix to his Submatrix(i,j).
                void Submatrix(unsigned int i, unsigned int j)
                {
                        MatrixTmpl<T> aux( GetDelRow(i) );
                        (*this) = aux.GetDelColumn(j);
                }

                // Solving linear equations [A][x]=[b]
                // friend MatrixTmpl<T> Solve(const MatrixTmpl<T>& A, const MatrixTmpl<T>& b);

                // LU Decomposition
                // Decompose the matrix into a product of lower and upper triangular matrices.
                // Decompose();
                // Compute forward/backward substitution on decomposed LU matrix product
                // Substitution(LU);

                // Eigenvalues and Eigen 
                // friend MatrixTmpl<T> Eigenvalues(const MatrixTmpl<T>& m);
                // friend MatrixTmpl<T> Eigenvectors(const MatrixTmpl<T>& m);

                // Gram Schmidt Orthonormalization
                // friend MatrixTmpl<T> Orth(const MatrixTmpl<T>& m);

                /**** Get/Set an element *****/
                void SetAt (unsigned int iPosition, unsigned int jPosition, T element)
                {
                        CLAM_ASSERT(iPosition < mNumRows,
                                "Index beyond matrix dimension");
                        CLAM_ASSERT(iPosition >= 0,
                                "Index beyond matrix dimension");
                        CLAM_ASSERT(jPosition < mNumColumns,
                                "Index beyond matrix dimension");
                        CLAM_ASSERT(jPosition >= 0,
                                "Index beyond matrix dimension");

                        MatrixBuffer()[mNumColumns*iPosition+jPosition] = element;
                }

                T GetAt (unsigned int iPosition, unsigned int jPosition) const
                {
                        CLAM_ASSERT(iPosition < mNumRows,
                                "Index beyond matrix dimension");
                        CLAM_ASSERT(iPosition >= 0,
                                "Index beyond matrix dimension");
                        CLAM_ASSERT(jPosition < mNumColumns,
                                "Index beyond matrix dimension");
                        CLAM_ASSERT(jPosition >= 0,
                                "Index beyond matrix dimension");

                        return MatrixBuffer()[mNumColumns*iPosition+jPosition];
                }

                // Get one column
                friend MatrixTmpl<T> GetColumn(unsigned int column, MatrixTmpl<T>& m )
                {
                        CLAM_ASSERT(column<m.mNumColumns,
                                "Column beyond matrix dimensions");
                        CLAM_ASSERT(column >= 0,
                                "Column beyond matrix dimensions");

                        MatrixTmpl<T> ret(m.mNumRows, 1); // Column vector
                        for (unsigned int i=0; i<m.mNumRows; i++)
                                ret(i,0) = m(i,column);
                        return ret;
                }

                // Get one row
                friend MatrixTmpl<T> GetRow(unsigned int row, MatrixTmpl<T>& m)
                {
                        CLAM_ASSERT(row < m.mNumRows,
                                "Row beyond matrix dimensions");
                        CLAM_ASSERT(row >= 0,
                                "Row beyond matrix dimensions");
        
                        MatrixTmpl<T> ret(1, m.mNumColumns); // Row vector
                        for (unsigned int i=0; i<m.mNumColumns; i++)
                                ret(0,i) = m(row,i);
                        return ret;
                }

                /* Apply an unary function to all the elements of the matrix */
                /* sin, cos, sqrt, cbrt, exp, log, log10, asin, acos, tan, atan */
                void Apply( T (*f)(T) )
                {
                        for (unsigned int i=0; i<mNumRows; i++)
                                for(unsigned int j=0; j<mNumColumns; j++)
                                        (*this)(i,j) = f( (*this)(i,j) );
                 }

                 void Apply( T (*f)(T,int),int parameter )
                {
                        for (unsigned int i=0; i<mNumRows; i++)
                                for(unsigned int j=0; j<mNumColumns; j++)
                                        (*this)(i,j) = f( (*this)(i,j), parameter );
                 }

                friend MatrixTmpl<T> GetApply( const MatrixTmpl<T> &m, double f(double) )
                {
                        MatrixTmpl<T> ret(m.mNumRows, m.mNumColumns);
                        for (unsigned int i=0; i<m.mNumRows; i++)
                                for(unsigned int j=0; j<m.mNumColumns; j++)
                                        ret(i,j) = T(f( m(i,j) ));
                        return ret;
                 }

                friend MatrixTmpl<T> AbsMatrix(const MatrixTmpl<T> &m) // absolute value
                {
                        MatrixTmpl<T> ret(m.mNumRows, m.mNumColumns);
                        for (unsigned int i=0; i<m.mNumRows; i++)
                                for(unsigned int j=0; j<m.mNumColumns; j++)
                                        ret(i,j) = abs( m(i,j) );
                        return ret;
                }

                /**** Operators ****/
                T& operator () (unsigned int iPosition,unsigned int jPosition) const {
                        CLAM_ASSERT(iPosition < mNumRows,
                                "Index beyond matrix dimension");
                        CLAM_ASSERT(iPosition >= 0,
                                "Index beyond matrix dimension");
                        CLAM_ASSERT(jPosition < mNumColumns,
                                "Index beyond matrix dimension");
                        CLAM_ASSERT(jPosition >= 0,
                                "Index beyond matrix dimension");

                        return ( MatrixBuffer()[mNumColumns*iPosition+jPosition] ); 
                }

                const MatrixTmpl<T>& operator = (const MatrixTmpl<T>& originalMatrix)
                {
                        // MatrixBuffer() = originalMatrix.MatrixBuffer();
                        *mpMatrixBuffer = *(originalMatrix.mpMatrixBuffer); 
                        //*mpMatrixBuffer = originalMatrix.MatrixBuffer();
                        mNumRows = originalMatrix.mNumRows;
                        mNumColumns = originalMatrix.mNumColumns;
                        return *this;
                }

                const MatrixTmpl<T>& operator = (const T element)
                {
                        MatrixBuffer().SetSize(1);
                        MatrixBuffer()[0] = element;
                        mNumRows = mNumColumns = 1;
                        return *this;
                }

                const MatrixTmpl<T>& operator += (const MatrixTmpl<T>& newMatrix)
                {
                        // Adding element by element
                        CLAM_ASSERT(mNumRows == newMatrix.mNumRows,"Adding matrix with different dimensions");
                        CLAM_ASSERT(mNumColumns == newMatrix.mNumColumns,"Adding matrix with different dimensions");
                        for (unsigned int i=0; i< (mNumRows*mNumColumns) ; i++)
                                MatrixBuffer()[i] += newMatrix.MatrixBuffer()[i];
                        return *this;
                }

                const MatrixTmpl<T>& operator -= (const MatrixTmpl<T>& newMatrix)
                {
                        // Substract element by element
                        CLAM_ASSERT(mNumRows == newMatrix.mNumRows,"Substracting matrix with different dimensions");
                        CLAM_ASSERT(mNumColumns == newMatrix.mNumColumns,"Substracting matrix with different dimensions");
                        for (unsigned int i=0; i< (mNumRows*mNumColumns) ; i++)
                                MatrixBuffer()[i] -= newMatrix.MatrixBuffer()[i];
                        return *this;
                }

                friend MatrixTmpl<T> operator + (MatrixTmpl<T>& m1, MatrixTmpl<T>& m2)
                {
                        CLAM_ASSERT(m1.mNumRows == m2.mNumRows,"Adding matrix with different dimensions");
                        CLAM_ASSERT(m1.mNumColumns == m2.mNumColumns,"Adding matrix with different dimensions");
                        MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns); // Construction of an Vector object
                        for (unsigned int i=0; i<ret.mNumRows; i++)
                                for (unsigned int j=0; j<ret.mNumColumns; j++)
                                        ret(i,j) = m1(i,j) + m2(i,j);
                        return ret;
                }

                // add an element to all the matrix elements
                friend MatrixTmpl<T> operator + (const MatrixTmpl<T>& m1, const T & element)
                {
                        MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns); 
                        for (unsigned int i=0; i<ret.mNumRows; i++)
                                for (unsigned int j=0; j<ret.mNumColumns; j++)
                                        ret(i,j) = m1(i,j) + element;
                        return ret;
                }

                friend MatrixTmpl<T> operator - (MatrixTmpl<T>& m1, MatrixTmpl<T>& m2)
                {
                        CLAM_ASSERT(m1.mNumRows == m2.mNumRows,"Substracting matrix with different dimensions");
                        CLAM_ASSERT(m1.mNumColumns == m2.mNumColumns,"Substracting matrix with different dimensions");
                        MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns); // Construction of an Vector object
                        for (unsigned int i=0; i<ret.mNumRows; i++)
                                for (unsigned int j=0; j<ret.mNumColumns; j++)
                                        ret(i,j) = m1(i,j) - m2(i,j);
                        return ret;
                }

                // substract an element to all the matrix elements
                friend MatrixTmpl<T> operator - (MatrixTmpl<T>& m1, const T element)
                {
                        MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns); 
                        for (unsigned int i=0; i<ret.mNumRows; i++)
                                for (unsigned int j=0; j<ret.mNumColumns; j++)
                                        ret(i,j) = m1(i,j) - element;
                        return ret;
                }

                friend MatrixTmpl<T> operator * (T scalar,const MatrixTmpl<T>& m)
                {
                        MatrixTmpl<T> mult(m.mNumRows, m.mNumColumns);
                        for (unsigned int i=0; i<(mult.mNumRows*mult.mNumColumns); i++)
                                mult.MatrixBuffer()[i] = scalar * m.MatrixBuffer()[i];
                        return mult;
                }

                friend MatrixTmpl<T> operator * (const MatrixTmpl<T>& m1, const MatrixTmpl<T>& m2)
                {
                        CLAM_ASSERT( m1.mNumColumns == m2.mNumRows,"Multiplications with incompatible arrays");
                        MatrixTmpl<T> ret(m1.mNumRows, m2.mNumColumns);
                        for (unsigned int i=0; i<ret.mNumRows; i++)
                                for (unsigned int j=0; j<ret.mNumColumns; j++) {
                                        ret(i,j) = 0;
                                        for( unsigned int k=0; k<m1.mNumColumns; k++)
                                                ret(i,j) += m1(i,k)*m2(k,j);
                                }
                        return ret;
                }

                friend MatrixTmpl<T> operator / (const MatrixTmpl<T>& m, T scalar)
                {
                        MatrixTmpl<T> ret(m.mNumRows, m.mNumColumns);
                        for (unsigned int i=0; i<(ret.mNumRows*ret.mNumColumns); i++)
                                ret.MatrixBuffer()[i] = m.MatrixBuffer()[i] / scalar;
                        return ret;
                }

                friend bool operator == (const MatrixTmpl<T>& m1, const MatrixTmpl<T>& m2)
                {
                        if ( (m1.mNumRows == m2.mNumRows) && (m1.mNumColumns == m2.mNumColumns) && (m1.MatrixBuffer() ==
        m2.MatrixBuffer()) )
                                return true;
                        else
                                return false; 
                }


                Array <T>& MatrixBuffer() const { return *mpMatrixBuffer;}

        protected:

                // Dimensions
                unsigned int mNumRows;
                unsigned int mNumColumns;

                // Buffer
                Array <T>* mpMatrixBuffer;

        };

        template <class T>
        std::istream& operator >> (std::istream & stream, MatrixTmpl<T> & a);
        
        template <class T>
        std::ostream& operator << (std::ostream & stream, const MatrixTmpl<T> & a);
        
} // namespace CLAM

#endif // _MatrixTmplDec_