HMSBEAGLE
1.0.0
|
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