|| || || || || || || ||

---

MatrixTem.h

Go to the documentation of this file.

#ifndef Impala_Core_Matrix_Matrix_h
#define Impala_Core_Matrix_Matrix_h

#include <iostream>
#include "Core/Array/Element/ArithTypes.h"
#include "Core/Vector/VectorTem.h"

namespace Impala
{
namespace Core
{
namespace Matrix
{

using namespace Impala::Core::Vector;


template<class C>
class MatrixTem {
public:


        MatrixTem()
        {
                _nr = 0;
                _nc = 0;
                _data = 0;
        }

    MatrixTem(int nRow, int nCol)
        {
                _nr = nRow;
                _nc = nCol;
                _data = new C[_nr * _nc];
        }

        MatrixTem(int nRow, int nCol, C a)
        {
                _nr = nRow;
                _nc = nCol;
                _data = new C[_nr * _nc];

                C* t = _data;
                int i = nElem();
                while (--i >= 0)
                        *t++ = a;
        }

        MatrixTem(int nRow, int nCol, C* data)
        {
                _nr = nRow;
                _nc = nCol;
                _data = data;
        }

        MatrixTem(const MatrixTem &m)
        {
                _nr = m.nRow();
                _nc = m.nCol();
                _data = new C[_nr * _nc];

                C* t = m._data;
                C* u = _data;
                int i = m.nElem();
                while (--i >= 0)
                        *u++ = *t++;
        }

        //MatrixTem(const VectorTem<C>& v)


        // Destructor
        ~MatrixTem()
        {
                //delete[] _data ;
        }


    int
        nRow() const
        {
                return _nr;
        }

    int
        nCol() const
        {
                return _nc;
        }

    int
        nElem() const
        {
                return _nr * _nc;
        }

    int
        valid() const
        {
                return ((_nc != 0) && (_nr != 0));
        }

        // DIRTY : testing only :-)
        C*
        dataPtr() const
        {
                return _data;
        }




        MatrixTem<C>&
        operator=(const MatrixTem<C>& m)
        {
                if (this != &m) {
                        delete [] _data;
                        _nr = m.nRow();
                        _nc = m.nCol();
                        _data = new C [_nr * _nc];
                        C *t = _data;
                        C *u = m._data;
                        int i = m.nElem();
                        while (--i >= 0)
                                *t++ = *u++;
                }
                return *this;
        }

        C*
        operator[](int i) const
        {
                return &_data[i*_nc];
        }

        MatrixTem<C>
        operator*(double          val) const
        {
                MatrixTem<C> m(*this);
                C* t = m._data;
                int i = nElem();
                while (--i >= 0)
                        *t++ = *t * val ;
                return m;
        }

        MatrixTem<C>
        operator*(const MatrixTem<C>& a) const
        {
                if (nCol() != a.nRow()) {
                        std::cerr << "nonconformant MatrixTem * MatrixTem operands."
                      << std::endl;
                        return MatrixTem<C>(0,0);
                }
                MatrixTem<C> m((*this).nRow(), a.nCol(), 0.0);
                double sum;
                int i, j, k;
                for (i=0 ; i<nRow() ; i++) {
                        for (j=0 ; j<a.nCol() ; j++) {
                                sum = 0;
                                for (k=0 ; k<nCol() ; k++)
                                        sum += (*this)[i][k] * a[k][j];
                                m[i][j] = sum;
                        }
                }
                return m;
        }

    template<class CC>
        friend VectorTem<CC>
        operator*(const VectorTem<CC>& , const MatrixTem<CC>& );

    template<class CC>
        friend VectorTem<CC>
        operator*(const MatrixTem<CC>& , const VectorTem<CC>& );

        MatrixTem<C>
        operator+(const MatrixTem<C>& a) const 
{
    if ((a.nCol() != nCol()) || (a.nRow() != nRow()) ) {
        std::cerr << "nonconformant MatrixTem + MatrixTem operands." << std::endl;
        return MatrixTem<C>(0,0);
    }
    MatrixTem<C> m(nRow(), nCol());
    int i, j;
    for (i=0 ; i<nRow() ; i++) {
        for (j=0 ; j<nCol() ; j++)
            m[i][j] = a[i][j] + (*this)[i][j];
    }
    return m;
}

        MatrixTem<C>            operator+(double          val) const 
{

    MatrixTem<C> m(*this);
    C* t = m._data;
    int i = nElem();
    while (--i >= 0)
        *t++ = val + *t;
    return m;
}



VectorTem<C>
getRow(long r) const
{

        if(r>= _nr || r<0){
        std::cerr << "getRow: rownr exceeded the number of rows in MatrixTem." 
                  << std::endl;
        return VectorTem<C>(0);
    }

        VectorTem<C> tmp(_nc);

        for(int i = 0; i<_nc; i++)
                tmp[i] = (*this)[r][i] ;

        return tmp ;
}

VectorTem<C>
getCol(long c) const
{

        if(c>= _nc || c<0){
        std::cerr << "getCol: Colnr exceeded number of columns in MatrixTem."
                  << std::endl;
        return VectorTem<C>(0);
    }

        VectorTem<C> tmp(_nr);
        for(int i = 0; i<_nr; i++)
        {
                tmp[i] = (*this)[i][c] ;
        }
        return tmp;

}


MatrixTem<C> 
setRow(long r, VectorTem<C> v, bool makeCopy = false) 
{

        if(r>= _nr || r<0){
        std::cerr << "setRow: rownr exceeded the number of rows in MatrixTem."
                  << std::endl;
        return MatrixTem<C>(0,0) ;
    }
        if(_nc != v.Size()){
        std::cerr << "setRow: v.nElem not equal to matrix dimension" << std::endl;
        return MatrixTem<C>(0,0) ;
    }

        MatrixTem<C> tmp(0,0);

        if(makeCopy) {

                tmp = MatrixTem<C>(_nr,_nc);

                for(int row = 0; row<_nr; row++)
                        for(int col = 0; col<_nc; col++)
                                if(row!=r) tmp[row][col] = (*this)[row][col] ;
                                else tmp[row][col] = v[col] ;
        } 
        else { // only change the original, faster

                for(int col = 0; col<_nc; col++)
                        (*this)[r][col] = v[col] ;
        }


        return tmp;

}

MatrixTem<C> 
setCol(long c, VectorTem<C> v, bool makeCopy = false) 
{

        if(c>= _nc || c<0){
        std::cerr << "setCol: Colnr exceeded number of columns in MatrixTem."
                  << std::endl;
        return MatrixTem<C>(0,0) ;
    }
        if(_nr != v.Size()){
        std::cerr << "setCol: v.nElem not equal to matrix dimension" << std::endl;
        return MatrixTem<C>(0,0) ;
    }

        MatrixTem<C> tmp(0,0);

        if(makeCopy) {

                tmp = MatrixTem<C>(_nr,_nc);

                for(int row = 0; row<_nr; row++)
                        for(int col = 0; col<_nc; col++)
                                if(col!=c) tmp[row][col] = (*this)[row][col] ;
                                else tmp[row][col] = v[row] ;
        }
        else { // only change the original, faster

                for(int row = 0; row<_nr; row++)
                        (*this)[row][c] = v[row] ;
        }

        return tmp;

}

    MatrixTem<C>                i() const;

        MatrixTem<C>
        t() const
        {
                MatrixTem<C> m(nCol(), nRow());
                int i, j;
                for (i=0 ; i<nRow() ; i++)
                        for (j=0 ; j<nCol() ; j++)
                                m[j][i] = (*this)[i][j];
                return m;
        }


    MatrixTem<C>
        covariance(VectorTem<C> &mean) const ;

    MatrixTem<C>
        svd(VectorTem<C>& W, MatrixTem<C>& V) const;

    MatrixTem<C>
        klm(int newDim) const;


//private:
        friend class VectorTem<C>;

    int         _nr;    // number of rows
    int         _nc;    // number of columns
    C*          _data;  // data pointer

};


//typedef MatrixTem<Real64>  Matrix;

template<class C>
inline VectorTem<C>
operator*(const VectorTem<C> &a, const MatrixTem<C> &b)
{
    if (a.Size() != b.nRow()) {
        std::cerr << "nonconformant VectorTem<C> * MatrixTem operands." << std::endl;
        return VectorTem<C>(0);
    }
    VectorTem<C> v(b.nCol());
    double sum;
    int i, j;
    for (i=0 ; i<b.nCol() ; i++) {
        sum = 0;
        for (j=0 ; j<b.nRow() ; j++)
            sum += a[j] * b[j][i];
        v[i] = sum;
    }
    return v;
}


template<class C>
inline VectorTem<C>
operator*(const MatrixTem<C> &a, const VectorTem<C> &b)
{
    if (b.Size() != a.nCol()) {
        std::cerr << "nonconformant MatrixTem * VectorTem<C> operands." << std::endl;
        return VectorTem<C>(0);
    }
    VectorTem<C> v(a.nRow());
    double sum;
    int i, j;
    for (i=0 ; i<a.nRow() ; i++) {
        sum = 0;
        for (j=0 ; j<a.nCol() ; j++)
            sum += a[i][j] * b[j];
        v[i] = sum;
    }
    return v;
}


template<class C> 
static void svdcmp(C* a, VectorTem<C> &w, C* v, int n, int m);

template<class C> 
static int ludcmp(C* a, int n, short *indx, double *d);

template<class C> 
static void lubksb(C* a, int n, short *indx, double *b);

template<class C> MatrixTem<C> 
MatrixTem<C>::i() const
{
    if (nRow() != nCol()) {
        std::cerr << "Inverse: matrix is not square!." << std::endl;
        return MatrixTem<C>(0, 0);
    }

    int size = nRow();
    MatrixTem<C> m(size,size), tmp(*this);
    short* idx = new short [size];
    double d;

    if (!ludcmp(tmp._data, size, idx, &d)) {
        std::cerr << "Inverse: singular matrix can't be inverted." << std::endl;
        delete [] idx;
        return MatrixTem<C>(0,0);
    }
    
    double * res = new double [size];

    for (int j = 0; j < size; j++) {
        for (int i = 0; i < size; i++)
            res[i] = 0.0;
        res[j] = 1.0;
        lubksb(tmp._data, size, idx, res);
        for (int i = 0; i < size; i++)
            m[i][j] = res[i];
    }

     delete [] res;
     delete [] idx;

     return m;
}

#ifndef HxMatrix_EPS_DEF 
#define HxMatrix_EPS_DEF
const double HxMatrix_EPS = 1e-12;
#endif

// ===================== matrix inverse =====================

template<class C> 
static int ludcmp(C* a, int n, short *indx, double *d)
{
    double *vv = new double [n];

    *d=1.0;
        int i, j, k;
    for ( i=0;i<n;i++) {
        double big=0.0;
        for (j=0;j<n;j++) {
            double temp = fabs(a[i*n+j]);
            if (temp > big)
                big=temp;
        }
        if (big == 0.0)
            return(0); // Singular matrix 
        vv[i]=1.0/big;
    }
    for (j=0;j<n;j++) {
        for (i=0;i<j;i++) {
            double sum=a[i*n+j];
            for (k=0;k<i;k++)
                sum -= a[i*n+k]*a[k*n+j];
            a[i*n+j]=sum;
        }
        double big=0.0;
        int imax = 0;
        for (i=j;i<n;i++) {
            double sum=a[i*n+j];
            for (k=0;k<j;k++)
                sum -= a[i*n+k]*a[k*n+j];
            a[i*n+j]=sum;

            double dum = vv[i]*fabs(sum);
            if ( dum >= big) {
                big=dum;
                imax=i;
            }
        }
        if (j != imax) {
            for (k=0;k<n;k++) {
                double dum=a[imax*n+k];
                a[imax*n+k]=a[j*n+k];
                a[j*n+k]=dum;
            }
            *d = -(*d);
            vv[imax]=vv[j];
        }
        indx[j]=imax;
        if (a[j*n+j] == 0.0)
            a[j*n+j]=HxMatrix_EPS;
        if (j != n-1) {
            double dum=1.0/a[j*n+j];
            for (i=j+1;i<n;i++)
                a[i*n+j] *= dum;
        }
    }
    delete [] vv;

    return(1);
}

template<class C> 
static void lubksb(C* a, int n, short *indx, double *b)
{
    int ii = -1;

        int i;
    for (i=0;i<n;i++) {
        int ip=indx[i];
        double sum=b[ip];
        b[ip]=b[i];
        if (ii>=0)
            for (int j=ii;j<=i-1;j++)
                sum -= a[i*n+j]*b[j];
        else if (sum) ii=i;
        b[i]=sum;
    }
    for (i=n-1;i>=0;i--) {
        double sum=b[i];
        for (int j=i+1;j<n;j++)
            sum -= a[i*n+j]*b[j];
        b[i]=sum/a[i*n+i];
    }
}




// ===================== covariance matrix =====================
/*
 *  covariance matrix
 *  Author(s):
 *  Minh A. Hoang (minh@science.uva.nl)
 *  Jan van Gemert (jvgemert@science.uva.nl)
 * Compute the covariance matrix of the observed data
 * Output: [nc x nc] covariance matrix and a vector of mean-values for each dimension.
 */

template<class C> MatrixTem<C> 
MatrixTem<C>::covariance(VectorTem<C> &mean) const 
{
//      printf("Calculating the covariance matrix...");
        int i,j,k;

        MatrixTem<C> COVxy(_nc,_nc, 0.0) ;

        mean = VectorTem<C>(_nc) ; //means
        for (i=0; i<_nc; i++) mean[i] = 0;

        for (k=0; k<_nr; k++)
                for (i=0; i<_nc; i++) {
                        mean[i] += (*this)[k][i];
                        for (j=i; j<_nc; j++)
                                COVxy[i][j] += ((*this)[k][i] * (*this)[k][j]) ;
                }

        for (i=0; i<_nc; i++) mean[i] /= _nr;

        for (i=0; i<_nc; i++)
                for (j=i; j<_nc; j++) {
                        COVxy[i][j] = COVxy[i][j]/_nr - mean[i]*mean[j];
                        COVxy[j][i] = COVxy[i][j];
                }

//      printf("...done\n");
        return COVxy;
}


// ===================== Karhunen Loeve mapping (PCA) =====================
/*
 *  Karhunen Loeve mapping (PCA)
 *  Author(s):
 *  Minh A. Hoang (minh@science.uva.nl)
 *  Jan van Gemert (jvgemert@science.uva.nl)
 *  Principal Componant Analysis to reduce matrix to newDim dimensions
 */

// KLM helping functions

/*
 Householder reduction of a real, symmetric matrix a[0..n-1][0..n-1]
 Output: a is replaced by the orthogonal matrix Q effecting the
 transformation. d[0..n-1] returns the diagonal elements of the
 tridiagonal matrix, and e[0..n-1] the off-diagonal elements, with e[0]=0;
 */
template<class C> 
void tred2(MatrixTem<C> &a, VectorTem<C> &d, VectorTem<C> &e)
{
        int i,j,k,l;
        double scale,hh,h,g,f;

        int n = d.Size() ;

        for (i=n-1; i>0; i--)
        {
                l = i-1;
                h = scale = 0.0;
                if (l>0)
                {
                        for (k=0; k<=l; k++) scale += fabs(a[i][k]);
                        if (scale == 0.0) e[i] = a[i][l];
                        else
                        {
                                for (k=0; k<=l; k++)
                                {
                                        a[i][k] /= scale; //use scale a's for transformation
                                        h += a[i][k]*a[i][k];
                                }
                                f = a[i][l];
                                g = (f >=0.0 ? -sqrt(h) : sqrt(h));
                                e[i] = scale*g;
                                h -= f*g;
                                a[i][l] = f-g;
                                f = 0.0;
                                for (j=0; j<=l; j++)
                                {
                                        a[j][i] = a[i][j]/h;
                                        g = 0.0;
                                        for (k=0; k<=j; k++) g += a[j][k]*a[i][k];
                                        for (k=j+1; k<=l; k++) g += a[k][j]*a[i][k];
                                        e[j] = g/h; //form element of p in temporarily unused element of e
                                        f += e[j]*a[i][j];
                                }
                                hh = f/(h+h);
                                for (j=0; j<=l; j++)
                                {
                                        f = a[i][j];
                                        e[j] = g = e[j]-hh*f;
                                        for (k=0; k<=j; k++) a[j][k] -= (f*e[k]+g*a[i][k]);
                                }
                        }
                }
                else e[i] = a[i][l];
                d[i] = h;
        }

        d[0] = 0.0;
        e[0] = 0.0;

        for (i=0; i<n; i++)
        {
                l = i-1;
                if (d[i])
                {
                        for (j=0; j<=l; j++)
                        {
                                g = 0.0;
                                for (k=0; k<=l; k++) g += a[i][k]*a[k][j];
                                for (k=0; k<=l; k++) a[k][j] -= g*a[k][i];
                        }
                }
                d[i] = a[i][i];
                a[i][i] = 1.0;
                for (j=0; j<=l; j++) a[j][i] = a[i][j] = 0.0;
        }

}

#define SIGN(a, b) ((b) < 0.0 ? -fabs(a): fabs(a))

/*
 QL algorithm with implicit shifts, to determin the eigenvalues and eigenvectors
 of a real, symetric, tridiagonal matrix, of of a real, symmetric matrix previously
 reduced by tred2
 Input: d[0..n-1] contains the diagonal elements of the tridiagonal matrix
 e[0..n-1] inputs the subdiagonal elements of the tridiagonal matrix with e[0] arbitrary
 Output: d returns the eigenvalues, e is destroyed.
 If the eigenvectors of a tridiagonal matrix are desired, the matrix z[0..n-1][0..n-1]
 is input as the identity matrix. If the eigenvectors of a matrix that has been reduced
 by tred2 are required, then z is input as the matrix output by tred2.
 In either case, the kth column of z returns the normalized eigenvector corresponding to d[k]
 */
template<class C> 
void tqli(VectorTem<C> &d, VectorTem<C> &e, MatrixTem<C> &z)
{

        int m,l,iter,i,k;
        double s,r,p,g,f,dd,c,b;

        int n = e.Size() ;

        for (i=1; i<n; i++) e[i-1]=e[i]; //renumber
        e[n-1] = 0.0;


        for (l=0; l<n; l++)
        {
                iter=0;
                do
                {
                        for (m=l; m<n-1; m++)
                        {
                                dd = fabs(d[m]) + fabs(d[m+1]);
                                if ( (fabs(e[m]) + dd) == dd) break;
                        }
                        if (m != l)
                        {
                                if (iter++ == 30) printf("Too many iterations in tqli\n");
                                g = (d[l+1]-d[l])/(2.0*e[l]);
                                r = sqrt(g*g + 1.0);//=pythag(g,1.0);
                                g = d[m] - d[l] + e[l]/(g+SIGN(r,g));
                                s = c = 1.0;
                                p = 0.0;
                                for (i=m-1; i>=l; i--)
                                {
                                        f = s*e[i];
                                        b = c*e[i];
                                        e[i+1] = (r=sqrt(f*f+g*g));//=(r=pythag(f,g));
                                        if (r == 0.0)
                                        {
                                                d[i+1] -= p;
                                                e[m] = 0.0;
                                                break;
                                        }
                                        s = f/r;
                                        c = g/r;
                                        g = d[i+1]-p;
                                        r = (d[i]-g)*s + 2.0*c*b;
                                        d[i+1] = g+(p=s*r);
                                        g = c*r-b;

                                        for (k=0; k<n; k++)
                                        {
                                                f = z[k][i+1];
                                                z[k][i+1] = s*z[k][i] + c*f;
                                                z[k][i] = c*z[k][i] - s*f;
                                        }
                                }
                                if (r==0.0 && i>=l) continue;
                                d[l] -= p;
                                e[l] = g;
                                e[m] = 0.0;
                        }
                } while (m != l);
        }

        //eigenvalues d[0..n-1], eigenvectors z[0..n-1][0..n-1], sort them now
        int j;
        for (i=0; i<n-1; i++)
        {
                p = d[k=i];
                for (j=i+1; j<n; j++) if (d[j]>=p) p = d[k=j];
                if (k!=i) { //insert
                        d[k] = d[i];
                        d[i] = p;
                        //exchange corresponding eigenvectors
                        for (j=0; j<n; j++)
                        {
                                p = z[j][i];
                                z[j][i] = z[j][k];
                                z[j][k] = p;
                        }
                }
        }
}


//  Karhunen Loeve mapping (PCA)


template<class C> MatrixTem<C> 
MatrixTem<C>::klm(int newDim) const
{
        printf("fix this function (MatrixTem::klm) before use\n");
        if (newDim<1 || newDim>_nc)
        {
                printf("klm: new dimensionality is invalid, dimsize is set as 1...\n");
                newDim = 1;
        }

        MatrixTem<C> newfvec(_nr,newDim);

/*      MatrixTem<C> a = this->covariance() ;

        VectorTem<C> d(_nc) ;

        VectorTem<C> e(_nc) ;

        tred2(a,d,e);

        tqli(d,e,a);

        // d[k] is an eigenvalue and a[][k] is a corresponding eigenvector

        for (int k=0; k<_nr; k++)
                for (int i=0; i<newDim; i++) {
                        double vl = 0.0;
                        for (int j=0; j<_nc; j++) vl += ((*this)[k][j]*a[j][i]);
                        newfvec[k][i] = vl;
                }*/

        return newfvec;
}

} // namespace Matrix
} // namespace Core
} // namespace Impala

#endif

---