matrix2d.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:     Carlos Oscar S. Sorzano (coss@cnb.csic.es)
 *
 * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 *
 * 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
 *
 *  All comments concerning this program package may be sent to the
 *  e-mail address 'xmipp@cnb.csic.es'
 ***************************************************************************/

#ifndef CORE_MATRIX2D_H_
#define CORE_MATRIX2D_H_

#include <complex>
#include "xmipp_random_mode.h"
#include "xmipp_macros.h"
#include "xmipp_error.h"
#include "xmipp_strings.h"

class FileName;

#ifdef XMIPP_MMAP
#include <sys/mman.h>
#endif

template<typename T>
class Matrix1D;

// Forward declarations
template<typename T>
class Matrix2D;

template<typename T>
void ludcmp(const Matrix2D<T>& A, Matrix2D<T>& LU, Matrix1D< int >& indx, T& d);

template<typename T>
void lubksb(const Matrix2D<T>& LU, Matrix1D< int >& indx, Matrix1D<T>& b);

template<typename T>
void svdcmp(const Matrix2D< T >& a,
            Matrix2D< double >& u,
            Matrix1D< double >& w,
            Matrix2D< double >& v);

void svbksb(Matrix2D< double >& u,
            Matrix1D< double >& w,
            Matrix2D< double >& v,
            Matrix1D< double >& b,
            Matrix1D< double >& x);


void cholesky(const Matrix2D<double> &M, Matrix2D<double> &L);

void schur(const Matrix2D<double> &M, Matrix2D<double> &O, Matrix2D<double> &T);


#define MATRIX2D_ARRAY(m) ((m).mdata)

#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m) \
    for (size_t i=0; i<(m).mdimy; i++) \
        for (size_t j=0; j<(m).mdimx; j++)

#define MAT_ELEM(m,i,j) ((m).mdata[((size_t)i)*(m).mdimx+((size_t)j)])

#define MAT_XSIZE(m) ((m).mdimx)

#define MAT_YSIZE(m) ((m).mdimy)

#define MAT_SIZE(m) ((m).mdim)

#define dMij(m, i, j)  MAT_ELEM(m, i, j)

#define dMn(m, n)  ((m).mdata[(n)])

#define M3x3_BY_V3x1(a, M, b) { \
        spduptmp0 = dMn(M, 0) * XX(b) + dMn(M, 1) * YY(b) + dMn(M, 2) * ZZ(b); \
        spduptmp1 = dMn(M, 3) * XX(b) + dMn(M, 4) * YY(b) + dMn(M, 5) * ZZ(b); \
        spduptmp2 = dMn(M, 6) * XX(b) + dMn(M, 7) * YY(b) + dMn(M, 8) * ZZ(b); \
        XX(a) = spduptmp0; YY(a) = spduptmp1; ZZ(a) = spduptmp2; }

#define M3x3_BY_M3x3(A, B, C) { \
        spduptmp0 = dMn(B,0) * dMn(C,0) + dMn(B,1) * dMn(C,3) + dMn(B,2) * dMn(C,6); \
        spduptmp1 = dMn(B,0) * dMn(C,1) + dMn(B,1) * dMn(C,4) + dMn(B,2) * dMn(C,7); \
        spduptmp2 = dMn(B,0) * dMn(C,2) + dMn(B,1) * dMn(C,5) + dMn(B,2) * dMn(C,8); \
        spduptmp3 = dMn(B,3) * dMn(C,0) + dMn(B,4) * dMn(C,3) + dMn(B,5) * dMn(C,6); \
        spduptmp4 = dMn(B,3) * dMn(C,1) + dMn(B,4) * dMn(C,4) + dMn(B,5) * dMn(C,7); \
        spduptmp5 = dMn(B,3) * dMn(C,2) + dMn(B,4) * dMn(C,5) + dMn(B,5) * dMn(C,8); \
        spduptmp6 = dMn(B,6) * dMn(C,0) + dMn(B,7) * dMn(C,3) + dMn(B,8) * dMn(C,6); \
        spduptmp7 = dMn(B,6) * dMn(C,1) + dMn(B,7) * dMn(C,4) + dMn(B,8) * dMn(C,7); \
        spduptmp8 = dMn(B,6) * dMn(C,2) + dMn(B,7) * dMn(C,5) + dMn(B,8) * dMn(C,8); \
        dMn(A, 0) = spduptmp0; \
        dMn(A, 1) = spduptmp1; \
        dMn(A, 2) = spduptmp2; \
        dMn(A, 3) = spduptmp3; \
        dMn(A, 4) = spduptmp4; \
        dMn(A, 5) = spduptmp5; \
        dMn(A, 6) = spduptmp6; \
        dMn(A, 7) = spduptmp7; \
        dMn(A, 8) = spduptmp8; }

#define M2x2_BY_V2x1(a, M, b) { \
        spduptmp0 = dMn(M, 0) * XX(b) + dMn(M, 1) * YY(b); \
        spduptmp1 = dMn(M, 2) * XX(b) + dMn(M, 3) * YY(b); \
        XX(a) = spduptmp0; \
        YY(a) = spduptmp1; }

#define M2x2_BY_CT(M2, M1, k) { \
        dMn(M2, 0) = dMn(M1, 0) * k; \
        dMn(M2, 1) = dMn(M1, 1) * k; \
        dMn(M2, 2) = dMn(M1, 2) * k; \
        dMn(M2, 3) = dMn(M1, 3) * k; }

#define M3x3_BY_CT(M2, M1, k) { \
        dMn(M2, 0) = dMn(M1, 0) * k; \
        dMn(M2, 1) = dMn(M1, 1) * k; \
        dMn(M2, 2) = dMn(M1, 2) * k; \
        dMn(M2, 3) = dMn(M1, 3) * k; \
        dMn(M2, 4) = dMn(M1, 4) * k; \
        dMn(M2, 5) = dMn(M1, 5) * k; \
        dMn(M2, 6) = dMn(M1, 6) * k; \
        dMn(M2, 7) = dMn(M1, 7) * k; \
        dMn(M2, 8) = dMn(M1, 8) * k; }

#define M4x4_BY_CT(M2, M1, k) { \
        dMn(M2, 0) = dMn(M1, 0) * k; \
        dMn(M2, 1) = dMn(M1, 1) * k; \
        dMn(M2, 2) = dMn(M1, 2) * k; \
        dMn(M2, 3) = dMn(M1, 3) * k; \
        dMn(M2, 4) = dMn(M1, 4) * k; \
        dMn(M2, 5) = dMn(M1, 5) * k; \
        dMn(M2, 6) = dMn(M1, 6) * k; \
        dMn(M2, 7) = dMn(M1, 7) * k; \
        dMn(M2, 8) = dMn(M1, 8) * k; \
  dMn(M2, 9) = dMn(M1, 9) * k; \
  dMn(M2,10) = dMn(M1,10) * k; \
  dMn(M2,11) = dMn(M1,11) * k; \
  dMn(M2,12) = dMn(M1,12) * k; \
  dMn(M2,13) = dMn(M1,13) * k; \
  dMn(M2,14) = dMn(M1,14) * k; \
  dMn(M2,15) = dMn(M1,15) * k;}

#define M2x2_INV(Ainv, A) { \
  spduptmp0 = dMn(A,0) * dMn(A,3) - dMn(A,1) * dMn(A,2);\
  if (spduptmp0==0.0) \
   REPORT_ERROR(ERR_NUMERICAL,"2x2 matrix is not invertible"); \
        spduptmp0 = 1.0 / spduptmp0; \
        dMn(Ainv, 0) =  dMn(A,3); \
        dMn(Ainv, 1) = -dMn(A,1); \
        dMn(Ainv, 2) = -dMn(A,2); \
        dMn(Ainv, 3) =  dMn(A,0); \
        M2x2_BY_CT(Ainv, Ainv, spduptmp0); }

#define M3x3_INV(Ainv, A) { \
        dMn(Ainv, 0) =   dMn(A,8)*dMn(A,4)-dMn(A,7)*dMn(A,5); \
        dMn(Ainv, 1) = -(dMn(A,8)*dMn(A,1)-dMn(A,7)*dMn(A,2)); \
        dMn(Ainv, 2) =   dMn(A,5)*dMn(A,1)-dMn(A,4)*dMn(A,2); \
        dMn(Ainv, 3) = -(dMn(A,8)*dMn(A,3)-dMn(A,6)*dMn(A,5)); \
        dMn(Ainv, 4) =   dMn(A,8)*dMn(A,0)-dMn(A,6)*dMn(A,2); \
        dMn(Ainv, 5) = -(dMn(A,5)*dMn(A,0)-dMn(A,3)*dMn(A,2)); \
        dMn(Ainv, 6) =   dMn(A,7)*dMn(A,3)-dMn(A,6)*dMn(A,4); \
        dMn(Ainv, 7) = -(dMn(A,7)*dMn(A,0)-dMn(A,6)*dMn(A,1)); \
        dMn(Ainv, 8) =   dMn(A,4)*dMn(A,0)-dMn(A,3)*dMn(A,1); \
        spduptmp0 = dMn(A,0)*dMn(Ainv,0)+dMn(A,3)*dMn(Ainv,1)+dMn(A,6)*dMn(Ainv,2); \
  if (spduptmp0==0.0) \
   REPORT_ERROR(ERR_NUMERICAL,"3x3 matrix is not invertible"); \
     spduptmp0 = 1.0 / spduptmp0; \
        M3x3_BY_CT(Ainv, Ainv, spduptmp0); }

#define M4x4_INV(Ainv, A) { \
        dMn(Ainv, 0) =   dMn(A,5)*(dMn(A,10)*dMn(A,15)-dMn(A,11)*dMn(A,14))+\
             dMn(A,6)*(dMn(A,11)*dMn(A,13)-dMn(A,9) *dMn(A,15))+\
             dMn(A,7)*(dMn(A,9) *dMn(A,14)-dMn(A,10)*dMn(A,13));\
  dMn(Ainv, 1) =   dMn(A,1)*(dMn(A,11)*dMn(A,14)-dMn(A,10)*dMn(A,15))+\
       dMn(A,2)*(dMn(A,9) *dMn(A,15)-dMn(A,11)*dMn(A,13))+\
       dMn(A,3)*(dMn(A,10)*dMn(A,13)-dMn(A,9) *dMn(A,14));\
  dMn(Ainv, 2) =   dMn(A,1)*(dMn(A,6) *dMn(A,15)-dMn(A,7) *dMn(A,14))+\
       dMn(A,2)*(dMn(A,7) *dMn(A,13)-dMn(A,5) *dMn(A,15))+\
       dMn(A,3)*(dMn(A,5) *dMn(A,14)-dMn(A,6) *dMn(A,13));\
  dMn(Ainv, 3) =   dMn(A,1)*(dMn(A,7) *dMn(A,10)-dMn(A,6) *dMn(A,11))+\
       dMn(A,2)*(dMn(A,5) *dMn(A,11)-dMn(A,7) *dMn(A,9))+\
       dMn(A,3)*(dMn(A,6) *dMn(A,9) -dMn(A,5) *dMn(A,10));\
  dMn(Ainv, 4) =   dMn(A,4)*(dMn(A,11)*dMn(A,14)-dMn(A,10)*dMn(A,15))+\
       dMn(A,6)*(dMn(A,8) *dMn(A,15)-dMn(A,11)*dMn(A,12))+\
       dMn(A,7)*(dMn(A,10)*dMn(A,12)-dMn(A,8) *dMn(A,14));\
  dMn(Ainv, 5) =   dMn(A,0)*(dMn(A,10)*dMn(A,15)-dMn(A,11)*dMn(A,14))+\
       dMn(A,2)*(dMn(A,11)*dMn(A,12)-dMn(A,8) *dMn(A,15))+\
       dMn(A,3)*(dMn(A,8) *dMn(A,14)-dMn(A,10)*dMn(A,12));\
  dMn(Ainv, 6) =   dMn(A,0)*(dMn(A,7) *dMn(A,14)-dMn(A,6) *dMn(A,15))+\
       dMn(A,2)*(dMn(A,4) *dMn(A,15)-dMn(A,7) *dMn(A,12))+\
       dMn(A,3)*(dMn(A,6) *dMn(A,12)-dMn(A,4) *dMn(A,14));\
  dMn(Ainv, 7) =   dMn(A,0)*(dMn(A,6) *dMn(A,11)-dMn(A,7) *dMn(A,10))+\
       dMn(A,2)*(dMn(A,7) *dMn(A,8) -dMn(A,4) *dMn(A,11))+\
       dMn(A,3)*(dMn(A,4) *dMn(A,10)-dMn(A,6) *dMn(A,8));\
  dMn(Ainv, 8) =   dMn(A,4)*(dMn(A,9) *dMn(A,15)-dMn(A,11)*dMn(A,13))+\
       dMn(A,5)*(dMn(A,11)*dMn(A,12)-dMn(A,8) *dMn(A,15))+\
       dMn(A,7)*(dMn(A,8) *dMn(A,13)-dMn(A,9) *dMn(A,12));\
  dMn(Ainv, 9) =   dMn(A,0)*(dMn(A,11)*dMn(A,13)-dMn(A,9) *dMn(A,15))+\
       dMn(A,1)*(dMn(A,8) *dMn(A,15)-dMn(A,11)*dMn(A,12))+\
       dMn(A,3)*(dMn(A,9) *dMn(A,12)-dMn(A,8) *dMn(A,13));\
  dMn(Ainv,10) =   dMn(A,0)*(dMn(A,5) *dMn(A,15)-dMn(A,7) *dMn(A,13))+\
       dMn(A,1)*(dMn(A,7) *dMn(A,12)-dMn(A,4) *dMn(A,15))+\
       dMn(A,3)*(dMn(A,4) *dMn(A,13)-dMn(A,5) *dMn(A,12));\
  dMn(Ainv,11) =   dMn(A,0)*(dMn(A,7) *dMn(A,9) -dMn(A,5) *dMn(A,11))+\
       dMn(A,1)*(dMn(A,4) *dMn(A,11)-dMn(A,7) *dMn(A,8))+\
       dMn(A,3)*(dMn(A,5) *dMn(A,8) -dMn(A,4) *dMn(A,9));\
  dMn(Ainv,12) =   dMn(A,4)*(dMn(A,10)*dMn(A,13)-dMn(A,9) *dMn(A,14))+\
       dMn(A,5)*(dMn(A,8) *dMn(A,14)-dMn(A,10)*dMn(A,12))+\
       dMn(A,6)*(dMn(A,9) *dMn(A,12)-dMn(A,8) *dMn(A,13));\
  dMn(Ainv,13) =   dMn(A,0)*(dMn(A,9) *dMn(A,14)-dMn(A,10)*dMn(A,13))+\
       dMn(A,1)*(dMn(A,10)*dMn(A,12)-dMn(A,8) *dMn(A,14))+\
       dMn(A,2)*(dMn(A,8) *dMn(A,13)-dMn(A,9) *dMn(A,12));\
  dMn(Ainv,14) =   dMn(A,0)*(dMn(A,6) *dMn(A,13)-dMn(A,5) *dMn(A,14))+\
       dMn(A,1)*(dMn(A,4) *dMn(A,14)-dMn(A,6) *dMn(A,12))+\
       dMn(A,2)*(dMn(A,5) *dMn(A,12)-dMn(A,4) *dMn(A,13));\
  dMn(Ainv,15) =   dMn(A,0)*(dMn(A,5) *dMn(A,10)-dMn(A,6) *dMn(A,9))+\
       dMn(A,1)*(dMn(A,6) *dMn(A,8) -dMn(A,4) *dMn(A,10))+\
       dMn(A,2)*(dMn(A,4) *dMn(A,9) -dMn(A,5) *dMn(A,8));\
        spduptmp0 = dMn(A,0)*(dMn(A,5)*(dMn(A,10)*dMn(A,15)-dMn(A,11)*dMn(A,14))\
                              +dMn(A,6)*(dMn(A,11)*dMn(A,13)-dMn(A,9) *dMn(A,15))\
                              +dMn(A,7)*(dMn(A,9) *dMn(A,14)-dMn(A,10)*dMn(A,13)))\
                    +dMn(A,1)*(dMn(A,4)*(dMn(A,11)*dMn(A,14)-dMn(A,10)*dMn(A,15))\
                              +dMn(A,6)*(dMn(A,8) *dMn(A,15)-dMn(A,11)*dMn(A,12))\
                              +dMn(A,7)*(dMn(A,10)*dMn(A,12)-dMn(A,8) *dMn(A,14)))\
        +dMn(A,2)*(dMn(A,4)*(dMn(A,9) *dMn(A,15)-dMn(A,11)*dMn(A,13))\
         +dMn(A,5)*(dMn(A,11)*dMn(A,12)-dMn(A,8) *dMn(A,15))\
         +dMn(A,7)*(dMn(A,8) *dMn(A,13)-dMn(A,9) *dMn(A,12)))\
        +dMn(A,3)*(dMn(A,4)*(dMn(A,10)*dMn(A,13)-dMn(A,9) *dMn(A,14))\
         +dMn(A,5)*(dMn(A,8) *dMn(A,14)-dMn(A,10)*dMn(A,12))\
         +dMn(A,6)*(dMn(A,9) *dMn(A,12)-dMn(A,8) *dMn(A,13))); \
  if (spduptmp0==0.0) \
   REPORT_ERROR(ERR_NUMERICAL,"4x4 matrix is not invertible"); \
  spduptmp0 = 1.0 / spduptmp0; \
        M4x4_BY_CT(Ainv, Ainv, spduptmp0); }

template<typename T>
class Matrix2D
{
public:
    // The array itself
    T* mdata;

    // Destroy data
    bool destroyData;

    // Mapped data
    bool mappedData;

    // File descriptor for mapped files
    int fdMap;

    // Mapped data original pointer
    char* mdataOriginal;

    // Number of elements in X
    size_t mdimx;

    // Number of elements in Y
    size_t mdimy;

    // Total number of elements
    size_t mdim;



    Matrix2D()
    {
        coreInit();
    }

    Matrix2D(const FileName &fnMappedMatrix, int Ydim, int Xdim, size_t offset=0)
    {
        coreInit(fnMappedMatrix,Ydim,Xdim,offset);
    }

    Matrix2D(int Ydim, int Xdim)
    {
        coreInit();
        initZeros(Ydim, Xdim);
    }

    Matrix2D(const Matrix2D<T>& v)
    {
        coreInit();
        *this = v;
    }

    ~Matrix2D()
    {
        coreDeallocate();
    }

    Matrix2D<T>& operator=(const Matrix2D<T>& op1);



    void clear()
    {
        coreDeallocate();
        coreInit();
    }

    inline void convertTo(float out[3][3]) const {
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                out[i][j] = (*this)(i, j);
            }
        }
    }

    void coreInit(const FileName &fn, int Ydim, int Xdim, size_t offset=0);

    void coreInit();

    void coreAllocate( int _mdimy, int _mdimx);

    void coreDeallocate();



    void resize(size_t Ydim, size_t Xdim, bool noCopy=false);

    template<typename T1>
    void resize(const Matrix2D<T1> &v)
    {
        if (mdimx != v.mdimx || mdimy != v.mdimy)
            resize(v.mdimy, v.mdimx);
    }

    inline void resizeNoCopy(int Ydim, int Xdim)
    {
        resize(Ydim, Xdim, true);
    }

    template<typename T1>
    inline void resizeNoCopy(const Matrix2D<T1> &v)
    {
        if (mdimx != v.mdimx || mdimy != v.mdimy)
            resize(v.mdimy, v.mdimx, true);
    }

    void mapToFile(const FileName &fn, int Ydim, int Xdim, size_t offset=0);

    void submatrix(int i0, int j0, int iF, int jF);

    template <typename T1>
    inline bool sameShape(const Matrix2D<T1>& op) const
    {
        return ((mdimx == op.mdimx) && (mdimy == op.mdimy));
    }

    inline size_t Xdim() const
    {
        return mdimx;
    }

    inline size_t Ydim() const
    {
        return mdimy;
    }



    inline void initConstant(T val)
    {
        for (size_t j = 0; j < mdim; j++)
            mdata[j] = val;
    }

    inline void initConstant(size_t Ydim, size_t Xdim, T val)
    {
        if (mdimx!=Xdim || mdimy!=Ydim)
            resizeNoCopy(Ydim, Xdim);
        initConstant(val);
    }

    inline void initZeros()
    {
        memset(mdata,0,mdimx*mdimy*sizeof(T));
    }

    inline void initZeros(size_t Ydim, size_t Xdim)
    {
        if (mdimx!=Xdim || mdimy!=Ydim)
            resizeNoCopy(Ydim, Xdim);
        memset(mdata,0,mdimx*mdimy*sizeof(T));
    }

    template <typename T1>
    void initZeros(const Matrix2D<T1>& op)
    {
        if (mdimx!=op.mdimx || mdimy!=op.mdimy)
            resizeNoCopy(op);
        memset(mdata,0,mdimx*mdimy*sizeof(T));
    }

    void initRandom(size_t Ydim, size_t Xdim, double op1, double op2, RandomMode mode = RND_UNIFORM);

    void initGaussian(int Ydim, int Xdim, double op1=0., double op2=1.);

    inline void initIdentity()
    {
        initIdentity(MAT_XSIZE(*this));
    }

    inline void initIdentity(int dim)
    {
        initZeros(dim, dim);
        for (int i = 0; i < dim; i++)
            MAT_ELEM(*this,i,i) = 1;
    }
    void initGaussian(int dim, double var);



    inline T& operator()(int i, int j) const
    {
        return MAT_ELEM((*this),i,j);
    }
    inline void setVal(T val,int y, int x)
    {
        MAT_ELEM((*this),y,x)=val;
    }
    inline T getVal( int y, int x) const
    {
        return MAT_ELEM((*this),y,x);
    }

    inline Matrix2D<T> operator*(T op1) const
    {
        Matrix2D<T> tmp(*this);
        for (size_t i=0; i < mdim; i++)
            tmp.mdata[i] = mdata[i] * op1;
        return tmp;
    }

    inline Matrix2D<T> operator/(T op1) const
    {
        Matrix2D<T> tmp(*this);
        for (size_t i=0; i < mdim; i++)
            tmp.mdata[i] = mdata[i] / op1;
        return tmp;
    }

    inline friend Matrix2D<T> operator*(T op1, const Matrix2D<T>& op2)
    {
        Matrix2D<T> tmp(op2);
        for (size_t i=0; i < op2.mdim; i++)
            tmp.mdata[i] = op1 * op2.mdata[i];
        return tmp;
    }

    inline void operator*=(T op1)
    {
        for (size_t i=0; i < mdim; i++)
            mdata[i] *= op1;
    }

    inline void operator/=(T op1)
    {
        for (size_t i=0; i < mdim; i++)
            mdata[i] /= op1;
    }

    Matrix1D<T> operator*(const Matrix1D<T>& op1) const;

    Matrix2D<T> operator*(const Matrix2D<T>& op1) const;

    Matrix2D<T> operator+(const Matrix2D<T>& op1) const;

    void operator+=(const Matrix2D<T>& op1) const;

    Matrix2D<T> operator-(const Matrix2D<T>& op1) const;

    void operator-=(const Matrix2D<T>& op1) const;
    bool equal(const Matrix2D<T>& op,
               double accuracy = XMIPP_EQUAL_ACCURACY) const;
    bool equalAbs(const Matrix2D<T>& op,
               double accuracy = XMIPP_EQUAL_ACCURACY) const;



    T computeMax() const;

    T computeMin() const;

    void computeMaxAndMin(T &maxValue, T &minValue) const;

    void rowSum(Matrix1D<T> &sum) const;

        void colSum(Matrix1D<T> &sum) const;

    void rowEnergySum(Matrix1D<T> &sum) const;

    T** adaptForNumericalRecipes() const;

    inline T* adaptForNumericalRecipes2() const
    {
        return mdata - 1 - mdimx;
    }

    void loadFromNumericalRecipes(T** m, int Ydim, int Xdim);

    inline void killAdaptationForNumericalRecipes(T** m) const
    {
        free_Tmatrix(m, 1, mdimy, 1, mdimx);
    }

    void killAdaptationForNumericalRecipes2(T** m) const
        {}

    void read(const FileName &fn);

    void write(const FileName &fn) const;
    friend std::ostream& operator<<(std::ostream& ostrm, const Matrix2D<T>& v)
    {
        if (v.Xdim() == 0 || v.Ydim() == 0)
            ostrm << "NULL matrix\n";
        else
        {
            ostrm << std::endl;
            double max_val = v.computeMax();
            int prec = bestPrecision(max_val, 10);

            for (size_t i = 0; i < v.Ydim(); i++)
            {
                for (size_t j = 0; j < v.Xdim(); j++)
                {
                    ostrm << floatToString((double) v(i, j), 10, prec) << ' ';
                }
                ostrm << std::endl;
            }
        }

        return ostrm;
    }

    void fromVector(const Matrix1D<T>& op1);

    void toVector(Matrix1D<T>& op1) const;

    void copyToVector(std::vector<T> &v)
    {
        v.assign(mdata, mdata+mdim);
    }
    inline void copyFromVector(std::vector<T> &v,int Xdim, int Ydim)
    {
        if (mdimx!=Xdim || mdimy!=Ydim)
            resizeNoCopy(Ydim, Xdim);
        copy( v.begin(), v.begin()+v.size(), mdata);
    }

    void getRow(size_t i, Matrix1D<T>& v) const;

    void getCol(size_t j, Matrix1D<T>& v) const;

    void setRow(size_t i, const Matrix1D<T>& v);

    void setCol(size_t j, const Matrix1D<T>& v);

    void computeRowMeans(Matrix1D<double> &Xmr) const;

    void computeColMeans(Matrix1D<double> &Xmr) const;

    inline void setConstantCol(size_t j, T v)
    {
        if (mdimx == 0 || mdimy == 0)
            REPORT_ERROR(ERR_MATRIX_EMPTY, "setCol: Target matrix is empty");

        if (j>= mdimx)
            REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS, "setCol: Matrix subscript (j) out of range");

        for (size_t i = 0; i < mdimy; i++)
            MAT_ELEM(*this,i, j) = v;
    }

    void getDiagonal(Matrix1D<T> &d) const;

    T trace() const;

    T det() const;

    inline T det3x3() const
    {
        return (
                   dMij(*this,0,0)*( dMij(*this,2,2)*dMij(*this,1,1)-
                                     dMij(*this,2,1)*dMij(*this,1,2) )-
                   dMij(*this,1,0)*( dMij(*this,2,2)*dMij(*this,0,1)-
                                     dMij(*this,2,1)*dMij(*this,0,2) )+
                   dMij(*this,2,0)*( dMij(*this,1,2)*dMij(*this,0,1)-
                                     dMij(*this,1,1)*dMij(*this,0,2) )
               );
    }

    inline double norm()
    {
        double sum=0.;
        FOR_ALL_ELEMENTS_IN_MATRIX2D(*this)
        {
            T aux=MAT_ELEM(*this,i,j);
            sum+=aux*aux;
        }
        return sqrt(sum);
    }

    Matrix2D<T> transpose() const;

    void inv(Matrix2D<T>& result) const;

    void invAlgLib(Matrix2D<T>& result, bool use_lu=false) const;

    void svd(Matrix2D<double> &U, Matrix1D<double> &W, Matrix2D<double> &V) const
    {
      svdcmp(*this, U, W, V);
    }

    void eigs(Matrix2D<double> &U, Matrix1D<double> &W, Matrix2D<double> &V, Matrix1D<int> &indexes) const;

    Matrix2D<T> inv() const
    {
        Matrix2D<T> result;
        inv(result);

        return result;
    }

    void selfInverse()
    {
        Matrix2D<T> auxMatrix(*this);
        auxMatrix.inv(*this);
    }

    bool isIdentity() const;
};

typedef Matrix2D<double> DMatrix;
typedef Matrix2D<int> IMatrix;

template<typename T>
bool operator==(const Matrix2D<T>& op1, const Matrix2D<T>& op2)
{
    return op1.equal(op2);
}

template<typename T>
void ludcmp(const Matrix2D<T>& A, Matrix2D<T>& LU, Matrix1D< int >& indx, T& d);


template<typename T>
void lubksb(const Matrix2D<T>& LU, Matrix1D< int >& indx, Matrix1D<T>& b);

void svbksb(Matrix2D< double >& u,
            Matrix1D< double >& w,
            Matrix2D< double >& v,
            Matrix1D< double >& b,
            Matrix1D< double >& x);

//#define VIA_NR
#define VIA_BILIB

template<typename T>
void svdcmp(const Matrix2D< T >& a,
            Matrix2D< double >& u,
            Matrix1D< double >& w,
            Matrix2D< double >& v);

void generalizedEigs(const Matrix2D<double> &A, const Matrix2D<double> &B, Matrix1D<double> &D, Matrix2D<double> &P);

void firstEigs(const Matrix2D<double> &A, size_t M, Matrix1D<double> &D, Matrix2D<double> &P, bool Pneeded=true);

void lastEigs(const Matrix2D<double> &A, size_t M, Matrix1D<double> &D, Matrix2D<double> &P);

void eigsBetween(const Matrix2D<double> &A, size_t I1, size_t I2, Matrix1D<double> &D, Matrix2D<double> &P);

void allEigs(const Matrix2D<double> &A, std::vector< std::complex<double> > &eigs);

void connectedComponentsOfUndirectedGraph(const Matrix2D<double> &G, Matrix1D<int> &component);

template<typename T1, typename T2>
void typeCast(const Matrix2D<T1>& v1,  Matrix2D<T2>& v2)
{
    if (v1.mdim == 0)
    {
        v2.clear();
        return;
    }

    if (v1.mdimx!=v2.mdimx || v1.mdimy!=v2.mdimy)
        v2.resizeNoCopy(v1);
    for (unsigned long int n = 0; n < v1.mdim; n++)
        v2.mdata[n] = static_cast< T2 > (v1.mdata[n]);
}

template<typename T1>
void typeCast(const Matrix2D<T1>& v1,  Matrix2D<T1>& v2)
{
    v2=v1;
}

void orthogonalizeColumnsGramSchmidt(Matrix2D<double> &M);

void normalizeColumns(Matrix2D<double> &A);

void normalizeColumnsBetween0and1(Matrix2D<double> &A);

void subtractColumnMeans(Matrix2D<double> &A);

void matrixOperation_AtA(const Matrix2D <double> &A, Matrix2D<double> &B);

void matrixOperation_AAt(const Matrix2D <double> &A, Matrix2D<double> &C);

void matrixOperation_AB(const Matrix2D <double> &A, const Matrix2D<double> &B, Matrix2D<double> &C);

void matrixOperation_Ax(const Matrix2D <double> &A, const Matrix1D<double> &x, Matrix1D<double> &y);

void matrixOperation_ABt(const Matrix2D <double> &A, const Matrix2D <double> &B, Matrix2D<double> &C);

void matrixOperation_AtB(const Matrix2D <double> &A, const Matrix2D<double> &B, Matrix2D<double> &C);

void matrixOperation_Atx(const Matrix2D <double> &A, const Matrix1D<double> &x, Matrix1D<double> &y);

void matrixOperation_AtBt(const Matrix2D <double> &A, const Matrix2D<double> &B, Matrix2D<double> &C);

void matrixOperation_XtAX_symmetric(const Matrix2D<double> &X, const Matrix2D<double> &A, Matrix2D<double> &B);

void matrixOperation_IplusA(Matrix2D<double> &A);

void matrixOperation_IminusA(Matrix2D<double> &A);

void eraseFirstColumn(Matrix2D<double> &A);

void keepColumns(Matrix2D<double> &A, int j0, int jF);



#endif /* MATRIX2D_H_ */