HMSBEAGLE
1.0.0
|
00001 /* 00002 * EigenDecompositionCube.cpp 00003 * 00004 * Created on: Sep 24, 2009 00005 * Author: msuchard 00006 */ 00007 #ifndef _EigenDecompositionCube_hpp_ 00008 #define _EigenDecompositionCube_hpp_ 00009 00010 #include "libhmsbeagle/CPU/EigenDecompositionCube.h" 00011 00012 00013 namespace beagle { 00014 namespace cpu { 00015 00016 #if defined (BEAGLE_IMPL_DEBUGGING_OUTPUT) && BEAGLE_IMPL_DEBUGGING_OUTPUT 00017 const bool DEBUGGING_OUTPUT = true; 00018 #else 00019 const bool DEBUGGING_OUTPUT = false; 00020 #endif 00021 00022 BEAGLE_CPU_EIGEN_TEMPLATE 00023 EigenDecompositionCube<BEAGLE_CPU_EIGEN_GENERIC>::EigenDecompositionCube(int decompositionCount, 00024 int stateCount, 00025 int categoryCount, 00026 long flags) 00027 : EigenDecomposition<BEAGLE_CPU_EIGEN_GENERIC>(decompositionCount, 00028 stateCount, 00029 categoryCount, 00030 flags) { 00031 gEigenValues = (REALTYPE**) malloc(sizeof(REALTYPE*) * kEigenDecompCount); 00032 if (gEigenValues == NULL) 00033 throw std::bad_alloc(); 00034 00035 gCMatrices = (REALTYPE**) malloc(sizeof(REALTYPE*) * kEigenDecompCount); 00036 if (gCMatrices == NULL) 00037 throw std::bad_alloc(); 00038 00039 for (int i = 0; i < kEigenDecompCount; i++) { 00040 gCMatrices[i] = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount * kStateCount * kStateCount); 00041 if (gCMatrices[i] == NULL) 00042 throw std::bad_alloc(); 00043 00044 gEigenValues[i] = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount); 00045 if (gEigenValues[i] == NULL) 00046 throw std::bad_alloc(); 00047 } 00048 00049 matrixTmp = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount); 00050 firstDerivTmp = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount); 00051 secondDerivTmp = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount); 00052 } 00053 00054 BEAGLE_CPU_EIGEN_TEMPLATE 00055 EigenDecompositionCube<BEAGLE_CPU_EIGEN_GENERIC>::~EigenDecompositionCube() { 00056 00057 for(int i=0; i<kEigenDecompCount; i++) { 00058 free(gCMatrices[i]); 00059 free(gEigenValues[i]); 00060 } 00061 free(gCMatrices); 00062 free(gEigenValues); 00063 free(matrixTmp); 00064 free(firstDerivTmp); 00065 free(secondDerivTmp); 00066 } 00067 00068 BEAGLE_CPU_EIGEN_TEMPLATE 00069 void EigenDecompositionCube<BEAGLE_CPU_EIGEN_GENERIC>::setEigenDecomposition(int eigenIndex, 00070 const double* inEigenVectors, 00071 const double* inInverseEigenVectors, 00072 const double* inEigenValues) { 00073 00074 if (kFlags & BEAGLE_FLAG_INVEVEC_STANDARD) { 00075 int l = 0; 00076 for (int i = 0; i < kStateCount; i++) { 00077 gEigenValues[eigenIndex][i] = inEigenValues[i]; 00078 for (int j = 0; j < kStateCount; j++) { 00079 for (int k = 0; k < kStateCount; k++) { 00080 gCMatrices[eigenIndex][l] = inEigenVectors[(i * kStateCount) + k] 00081 * inInverseEigenVectors[(k * kStateCount) + j]; 00082 l++; 00083 } 00084 } 00085 } 00086 } else { 00087 int l = 0; 00088 for (int i = 0; i < kStateCount; i++) { 00089 gEigenValues[eigenIndex][i] = inEigenValues[i]; 00090 for (int j = 0; j < kStateCount; j++) { 00091 for (int k = 0; k < kStateCount; k++) { 00092 gCMatrices[eigenIndex][l] = inEigenVectors[(i * kStateCount) + k] 00093 * inInverseEigenVectors[k + (j*kStateCount)]; 00094 l++; 00095 } 00096 } 00097 } 00098 } 00099 00100 } 00101 00102 #define UNROLL 00103 00104 BEAGLE_CPU_EIGEN_TEMPLATE 00105 void EigenDecompositionCube<BEAGLE_CPU_EIGEN_GENERIC>::updateTransitionMatrices(int eigenIndex, 00106 const int* probabilityIndices, 00107 const int* firstDerivativeIndices, 00108 const int* secondDerivativeIndices, 00109 const double* edgeLengths, 00110 const double* categoryRates, 00111 REALTYPE** transitionMatrices, 00112 int count) { 00113 #ifdef UNROLL 00114 int stateCountModFour = (kStateCount / 4) * 4; 00115 #endif 00116 00117 if (firstDerivativeIndices == NULL && secondDerivativeIndices == NULL) { 00118 for (int u = 0; u < count; u++) { 00119 REALTYPE* transitionMat = transitionMatrices[probabilityIndices[u]]; 00120 int n = 0; 00121 for (int l = 0; l < kCategoryCount; l++) { 00122 00123 for (int i = 0; i < kStateCount; i++) { 00124 matrixTmp[i] = exp(gEigenValues[eigenIndex][i] * ((REALTYPE)edgeLengths[u] * categoryRates[l])); 00125 } 00126 00127 REALTYPE* tmpCMatrices = gCMatrices[eigenIndex]; 00128 for (int i = 0; i < kStateCount; i++) { 00129 for (int j = 0; j < kStateCount; j++) { 00130 REALTYPE sum = 0.0; 00131 #ifdef UNROLL 00132 int k = 0; 00133 for (; k < stateCountModFour; k += 4) { 00134 sum += tmpCMatrices[k + 0] * matrixTmp[k + 0]; 00135 sum += tmpCMatrices[k + 1] * matrixTmp[k + 1]; 00136 sum += tmpCMatrices[k + 2] * matrixTmp[k + 2]; 00137 sum += tmpCMatrices[k + 3] * matrixTmp[k + 3]; 00138 } 00139 for (; k < kStateCount; k++) { 00140 sum += tmpCMatrices[k] * matrixTmp[k]; 00141 } 00142 tmpCMatrices += kStateCount; 00143 #else 00144 for (int k = 0; k < kStateCount; k++) { 00145 sum += *tmpCMatrices++ * matrixTmp[k]; 00146 } 00147 #endif 00148 if (sum > 0) 00149 transitionMat[n] = sum; 00150 else 00151 transitionMat[n] = 0; 00152 n++; 00153 } 00154 if (T_PAD != 0) { 00155 transitionMat[n] = 1.0; 00156 n += T_PAD; 00157 } 00158 } 00159 } 00160 00161 if (DEBUGGING_OUTPUT) { 00162 int kMatrixSize = kStateCount * kStateCount; 00163 fprintf(stderr,"transitionMat index=%d brlen=%.5f\n", probabilityIndices[u], edgeLengths[u]); 00164 for ( int w = 0; w < (20 > kMatrixSize ? 20 : kMatrixSize); ++w) 00165 fprintf(stderr,"transitionMat[%d] = %.5f\n", w, transitionMat[w]); 00166 } 00167 } 00168 00169 00170 } else if (secondDerivativeIndices == NULL) { 00171 for (int u = 0; u < count; u++) { 00172 REALTYPE* transitionMat = transitionMatrices[probabilityIndices[u]]; 00173 REALTYPE* firstDerivMat = transitionMatrices[firstDerivativeIndices[u]]; 00174 int n = 0; 00175 for (int l = 0; l < kCategoryCount; l++) { 00176 00177 for (int i = 0; i < kStateCount; i++) { 00178 REALTYPE scaledEigenValue = gEigenValues[eigenIndex][i] * ((REALTYPE)categoryRates[l]); 00179 matrixTmp[i] = exp(scaledEigenValue * ((REALTYPE)edgeLengths[u])); 00180 firstDerivTmp[i] = scaledEigenValue * matrixTmp[i]; 00181 } 00182 00183 int m = 0; 00184 for (int i = 0; i < kStateCount; i++) { 00185 for (int j = 0; j < kStateCount; j++) { 00186 REALTYPE sum = 0.0; 00187 REALTYPE sumD1 = 0.0; 00188 for (int k = 0; k < kStateCount; k++) { 00189 sum += gCMatrices[eigenIndex][m] * matrixTmp[k]; 00190 sumD1 += gCMatrices[eigenIndex][m] * firstDerivTmp[k]; 00191 m++; 00192 } 00193 if (sum > 0) 00194 transitionMat[n] = sum; 00195 else 00196 transitionMat[n] = 0; 00197 firstDerivMat[n] = sumD1; 00198 n++; 00199 } 00200 if (T_PAD != 0) { 00201 transitionMat[n] = 1.0; 00202 firstDerivMat[n] = 0.0; 00203 n += T_PAD; 00204 } 00205 } 00206 } 00207 } 00208 } else { 00209 for (int u = 0; u < count; u++) { 00210 REALTYPE* transitionMat = transitionMatrices[probabilityIndices[u]]; 00211 REALTYPE* firstDerivMat = transitionMatrices[firstDerivativeIndices[u]]; 00212 REALTYPE* secondDerivMat = transitionMatrices[secondDerivativeIndices[u]]; 00213 int n = 0; 00214 for (int l = 0; l < kCategoryCount; l++) { 00215 00216 for (int i = 0; i < kStateCount; i++) { 00217 REALTYPE scaledEigenValue = gEigenValues[eigenIndex][i] * ((REALTYPE)categoryRates[l]); 00218 matrixTmp[i] = exp(scaledEigenValue * ((REALTYPE)edgeLengths[u])); 00219 firstDerivTmp[i] = scaledEigenValue * matrixTmp[i]; 00220 secondDerivTmp[i] = scaledEigenValue * firstDerivTmp[i]; 00221 } 00222 00223 int m = 0; 00224 for (int i = 0; i < kStateCount; i++) { 00225 for (int j = 0; j < kStateCount; j++) { 00226 REALTYPE sum = 0.0; 00227 REALTYPE sumD1 = 0.0; 00228 REALTYPE sumD2 = 0.0; 00229 for (int k = 0; k < kStateCount; k++) { 00230 sum += gCMatrices[eigenIndex][m] * matrixTmp[k]; 00231 sumD1 += gCMatrices[eigenIndex][m] * firstDerivTmp[k]; 00232 sumD2 += gCMatrices[eigenIndex][m] * secondDerivTmp[k]; 00233 m++; 00234 } 00235 if (sum > 0) 00236 transitionMat[n] = sum; 00237 else 00238 transitionMat[n] = 0; 00239 firstDerivMat[n] = sumD1; 00240 secondDerivMat[n] = sumD2; 00241 n++; 00242 } 00243 if (T_PAD != 0) { 00244 transitionMat[n] = 1.0; 00245 firstDerivMat[n] = 0.0; 00246 secondDerivMat[n] = 0.0; 00247 n += T_PAD; 00248 } 00249 } 00250 } 00251 } 00252 } 00253 } 00254 00255 } // cpu 00256 } // beagle 00257 00258 #endif // _EigenDecompositionCube_hpp_ 00259