HMSBEAGLE  1.0.0
libhmsbeagle/CPU/EigenDecompositionSquare.hpp
00001 /*
00002  * EigenDecompositionSquare.cpp
00003  *
00004  *  Created on: Sep 24, 2009
00005  *      Author: msuchard
00006  */
00007 #ifndef _EigenDecompositionSquare_hpp_
00008 #define _EigenDecompositionSquare_hpp_
00009 #include "EigenDecompositionSquare.h"
00010 #include "libhmsbeagle/beagle.h"
00011 
00012 //#if defined (BEAGLE_IMPL_DEBUGGING_OUTPUT) && BEAGLE_IMPL_DEBUGGING_OUTPUT
00013 //const bool DEBUGGING_OUTPUT = true;
00014 //#else
00015 //const bool DEBUGGING_OUTPUT = false;
00016 //#endif
00017 
00018 //#define DEBUG_COMPLEX
00019 
00020 namespace beagle {
00021 namespace cpu {
00022 
00023 BEAGLE_CPU_EIGEN_TEMPLATE
00024 EigenDecompositionSquare<BEAGLE_CPU_EIGEN_GENERIC>::EigenDecompositionSquare(int decompositionCount,
00025                                                                                                int stateCount,
00026                                                                                                int categoryCount,
00027                                                                                                long flags)
00028         : EigenDecomposition<BEAGLE_CPU_EIGEN_GENERIC>(decompositionCount,stateCount,categoryCount, flags) {
00029 
00030         isComplex = kFlags & BEAGLE_FLAG_EIGEN_COMPLEX;
00031 
00032         if (isComplex)
00033                 kEigenValuesSize = 2 * kStateCount;
00034         else
00035                 kEigenValuesSize = kStateCount;
00036 
00037     this->gEigenValues = (REALTYPE**) malloc(sizeof(REALTYPE*) * kEigenDecompCount);
00038     if (gEigenValues == NULL)
00039         throw std::bad_alloc();
00040 
00041     gEMatrices = (REALTYPE**) malloc(sizeof(REALTYPE*) * kEigenDecompCount);
00042     if (gEMatrices == NULL)
00043         throw std::bad_alloc();
00044 
00045     gIMatrices = (REALTYPE**) malloc(sizeof(REALTYPE*) * kEigenDecompCount);
00046        if (gIMatrices == NULL)
00047         throw std::bad_alloc();
00048 
00049     for (int i = 0; i < kEigenDecompCount; i++) {
00050         gEMatrices[i] = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount * kStateCount);
00051         if (gEMatrices[i] == NULL)
00052                 throw std::bad_alloc();
00053 
00054         gIMatrices[i] = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount * kStateCount);
00055         if (gIMatrices[i] == NULL)
00056                 throw std::bad_alloc();
00057 
00058         gEigenValues[i] = (REALTYPE*) malloc(sizeof(REALTYPE) * kEigenValuesSize);
00059         if (gEigenValues[i] == NULL)
00060                 throw std::bad_alloc();
00061     }
00062 
00063     matrixTmp = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount * kStateCount);
00064 }
00065 
00066 BEAGLE_CPU_EIGEN_TEMPLATE
00067 EigenDecompositionSquare<BEAGLE_CPU_EIGEN_GENERIC>::~EigenDecompositionSquare() {
00068 
00069         for(int i=0; i<kEigenDecompCount; i++) {
00070                 free(gEMatrices[i]);
00071                 free(gIMatrices[i]);
00072                 free(gEigenValues[i]);
00073         }
00074         free(gEMatrices);
00075         free(gIMatrices);
00076         free(gEigenValues);
00077         free(matrixTmp);
00078 }
00079     
00083 template<typename REALTYPE>
00084 void transposeSquareMatrix(REALTYPE* mat,
00085                            int size) {
00086     for (int i = 0; i < size - 1; i++) {
00087         for (int j = i + 1; j < size; j++) {
00088             REALTYPE tmp = mat[i * size + j];
00089             mat[i * size + j] = mat[j * size + i];
00090             mat[j * size + i] = tmp;
00091         }
00092     }
00093 }
00094 
00095 BEAGLE_CPU_EIGEN_TEMPLATE
00096 void EigenDecompositionSquare<BEAGLE_CPU_EIGEN_GENERIC>::setEigenDecomposition(int eigenIndex,
00097                                                                                              const double* inEigenVectors,
00098                                                      const double* inInverseEigenVectors,
00099                                                      const double* inEigenValues) {
00100     
00101         beagleMemCpy(gEigenValues[eigenIndex],inEigenValues,kEigenValuesSize);
00102         const int len = kStateCount * kStateCount;
00103         beagleMemCpy(gEMatrices[eigenIndex],inEigenVectors,len);
00104         beagleMemCpy(gIMatrices[eigenIndex],inInverseEigenVectors,len);
00105     if (kFlags & BEAGLE_FLAG_INVEVEC_TRANSPOSED) // TODO: optimize, might not need to transpose here
00106         transposeSquareMatrix(gIMatrices[eigenIndex], kStateCount);
00107 }
00108 
00109 BEAGLE_CPU_EIGEN_TEMPLATE
00110 void EigenDecompositionSquare<BEAGLE_CPU_EIGEN_GENERIC>::updateTransitionMatrices(int eigenIndex,
00111                                                         const int* probabilityIndices,
00112                                                         const int* firstDerivativeIndices,
00113                                                         const int* secondDerivativeIndices,
00114                                                         const double* edgeLengths,
00115                                                         const double* categoryRates,
00116                                                         REALTYPE** transitionMatrices,
00117                                                         int count) {
00118 
00119         const REALTYPE* Ievc = gIMatrices[eigenIndex];
00120         const REALTYPE* Evec = gEMatrices[eigenIndex];
00121         const REALTYPE* Eval = gEigenValues[eigenIndex];
00122         const REALTYPE* EvalImag = Eval + kStateCount;
00123     for (int u = 0; u < count; u++) {
00124         REALTYPE* transitionMat = transitionMatrices[probabilityIndices[u]];
00125         const double edgeLength = edgeLengths[u];
00126         int n = 0;
00127         for (int l = 0; l < kCategoryCount; l++) {
00128                         const REALTYPE distance = categoryRates[l] * edgeLength;
00129                 for(int i=0; i<kStateCount; i++) {
00130                         if (!isComplex || EvalImag[i] == 0) {
00131                                 const REALTYPE tmp = exp(Eval[i] * distance);
00132                                 for(int j=0; j<kStateCount; j++) {
00133                                         matrixTmp[i*kStateCount+j] = Ievc[i*kStateCount+j] * tmp;
00134                                 }
00135                         } else {
00136                                 // 2 x 2 conjugate block
00137                                 int i2 = i + 1;
00138                                 const REALTYPE b = EvalImag[i];
00139                                 const REALTYPE expat = exp(Eval[i] * distance);
00140                                 const REALTYPE expatcosbt = expat * cos(b * distance);
00141                                 const REALTYPE expatsinbt = expat * sin(b * distance);
00142                                 for(int j=0; j<kStateCount; j++) {
00143                                         matrixTmp[ i*kStateCount+j] = expatcosbt * Ievc[ i*kStateCount+j] +
00144                                                                               expatsinbt * Ievc[i2*kStateCount+j];
00145                                         matrixTmp[i2*kStateCount+j] = expatcosbt * Ievc[i2*kStateCount+j] -
00146                                                                                                       expatsinbt * Ievc[ i*kStateCount+j];
00147                                 }
00148                                 i++; // processed two conjugate rows
00149                         }
00150                 }
00151 
00152 #ifdef DEBUG_COMPLEX
00153                 fprintf(stderr,"[");
00154                 for(int i=0; i<16; i++)
00155                         fprintf(stderr," %7.5e,",matrixTmp[i]);
00156                 fprintf(stderr,"] -- complex debug\n");
00157                 exit(0);
00158 #endif
00159 
00160 
00161             for (int i = 0; i < kStateCount; i++) {
00162                 for (int j = 0; j < kStateCount; j++) {
00163                     REALTYPE sum = 0.0;
00164                     for (int k = 0; k < kStateCount; k++)
00165                         sum += Evec[i*kStateCount+k] * matrixTmp[k*kStateCount+j];
00166                     if (sum > 0)
00167                         transitionMat[n] = sum;
00168                     else
00169                         transitionMat[n] = 0;
00170                     n++;
00171                 }
00172 if (T_PAD != 0) {
00173                 transitionMat[n] = 1.0;
00174                 n += T_PAD;
00175 }
00176             }
00177         }
00178 
00179         if (DEBUGGING_OUTPUT) {
00180                 int kMatrixSize = kStateCount * kStateCount;
00181             fprintf(stderr,"transitionMat index=%d brlen=%.5f\n", probabilityIndices[u], edgeLengths[u]);
00182             for ( int w = 0; w < (20 > kMatrixSize ? 20 : kMatrixSize); ++w)
00183                 fprintf(stderr,"transitionMat[%d] = %.5f\n", w, transitionMat[w]);
00184         }
00185     }
00186 }
00187 
00188 }
00189 }
00190 
00191 #endif 
00192