7 #ifndef _EigenDecompositionCube_hpp_
8 #define _EigenDecompositionCube_hpp_
10 #include "libhmsbeagle/CPU/EigenDecompositionCube.h"
16 #if defined (BEAGLE_IMPL_DEBUGGING_OUTPUT) && BEAGLE_IMPL_DEBUGGING_OUTPUT
17 const bool DEBUGGING_OUTPUT =
true;
19 const bool DEBUGGING_OUTPUT =
false;
22 BEAGLE_CPU_EIGEN_TEMPLATE
23 EigenDecompositionCube<BEAGLE_CPU_EIGEN_GENERIC>::EigenDecompositionCube(
int decompositionCount,
27 : EigenDecomposition<BEAGLE_CPU_EIGEN_GENERIC>(decompositionCount,
31 gEigenValues = (REALTYPE**) malloc(
sizeof(REALTYPE*) * kEigenDecompCount);
32 if (gEigenValues == NULL)
33 throw std::bad_alloc();
35 gCMatrices = (REALTYPE**) malloc(
sizeof(REALTYPE*) * kEigenDecompCount);
36 if (gCMatrices == NULL)
37 throw std::bad_alloc();
39 for (
int i = 0; i < kEigenDecompCount; i++) {
40 gCMatrices[i] = (REALTYPE*) malloc(
sizeof(REALTYPE) * kStateCount * kStateCount * kStateCount);
41 if (gCMatrices[i] == NULL)
42 throw std::bad_alloc();
44 gEigenValues[i] = (REALTYPE*) malloc(
sizeof(REALTYPE) * kStateCount);
45 if (gEigenValues[i] == NULL)
46 throw std::bad_alloc();
49 matrixTmp = (REALTYPE*) malloc(
sizeof(REALTYPE) * kStateCount);
50 firstDerivTmp = (REALTYPE*) malloc(
sizeof(REALTYPE) * kStateCount);
51 secondDerivTmp = (REALTYPE*) malloc(
sizeof(REALTYPE) * kStateCount);
54 BEAGLE_CPU_EIGEN_TEMPLATE
55 EigenDecompositionCube<BEAGLE_CPU_EIGEN_GENERIC>::~EigenDecompositionCube() {
57 for(
int i=0; i<kEigenDecompCount; i++) {
59 free(gEigenValues[i]);
68 BEAGLE_CPU_EIGEN_TEMPLATE
69 void EigenDecompositionCube<BEAGLE_CPU_EIGEN_GENERIC>::setEigenDecomposition(
int eigenIndex,
70 const double* inEigenVectors,
71 const double* inInverseEigenVectors,
72 const double* inEigenValues) {
74 if (kFlags & BEAGLE_FLAG_INVEVEC_STANDARD) {
76 for (
int i = 0; i < kStateCount; i++) {
77 gEigenValues[eigenIndex][i] = inEigenValues[i];
78 for (
int j = 0; j < kStateCount; j++) {
79 for (
int k = 0; k < kStateCount; k++) {
80 gCMatrices[eigenIndex][l] = inEigenVectors[(i * kStateCount) + k]
81 * inInverseEigenVectors[(k * kStateCount) + j];
88 for (
int i = 0; i < kStateCount; i++) {
89 gEigenValues[eigenIndex][i] = inEigenValues[i];
90 for (
int j = 0; j < kStateCount; j++) {
91 for (
int k = 0; k < kStateCount; k++) {
92 gCMatrices[eigenIndex][l] = inEigenVectors[(i * kStateCount) + k]
93 * inInverseEigenVectors[k + (j*kStateCount)];
104 BEAGLE_CPU_EIGEN_TEMPLATE
105 void EigenDecompositionCube<BEAGLE_CPU_EIGEN_GENERIC>::updateTransitionMatrices(
int eigenIndex,
106 const int* probabilityIndices,
107 const int* firstDerivativeIndices,
108 const int* secondDerivativeIndices,
109 const double* edgeLengths,
110 const double* categoryRates,
111 REALTYPE** transitionMatrices,
114 int stateCountModFour = (kStateCount / 4) * 4;
117 if (firstDerivativeIndices == NULL && secondDerivativeIndices == NULL) {
118 for (
int u = 0; u < count; u++) {
119 REALTYPE* transitionMat = transitionMatrices[probabilityIndices[u]];
121 for (
int l = 0; l < kCategoryCount; l++) {
123 for (
int i = 0; i < kStateCount; i++) {
124 matrixTmp[i] = exp(gEigenValues[eigenIndex][i] * ((REALTYPE)edgeLengths[u] * categoryRates[l]));
127 REALTYPE* tmpCMatrices = gCMatrices[eigenIndex];
128 for (
int i = 0; i < kStateCount; i++) {
129 for (
int j = 0; j < kStateCount; j++) {
133 for (; k < stateCountModFour; k += 4) {
134 sum += tmpCMatrices[k + 0] * matrixTmp[k + 0];
135 sum += tmpCMatrices[k + 1] * matrixTmp[k + 1];
136 sum += tmpCMatrices[k + 2] * matrixTmp[k + 2];
137 sum += tmpCMatrices[k + 3] * matrixTmp[k + 3];
139 for (; k < kStateCount; k++) {
140 sum += tmpCMatrices[k] * matrixTmp[k];
142 tmpCMatrices += kStateCount;
144 for (
int k = 0; k < kStateCount; k++) {
145 sum += *tmpCMatrices++ * matrixTmp[k];
149 transitionMat[n] = sum;
151 transitionMat[n] = 0;
155 transitionMat[n] = 1.0;
161 if (DEBUGGING_OUTPUT) {
162 int kMatrixSize = kStateCount * kStateCount;
163 fprintf(stderr,
"transitionMat index=%d brlen=%.5f\n", probabilityIndices[u], edgeLengths[u]);
164 for (
int w = 0; w < (20 > kMatrixSize ? 20 : kMatrixSize); ++w)
165 fprintf(stderr,
"transitionMat[%d] = %.5f\n", w, transitionMat[w]);
170 }
else if (secondDerivativeIndices == NULL) {
171 for (
int u = 0; u < count; u++) {
172 REALTYPE* transitionMat = transitionMatrices[probabilityIndices[u]];
173 REALTYPE* firstDerivMat = transitionMatrices[firstDerivativeIndices[u]];
175 for (
int l = 0; l < kCategoryCount; l++) {
177 for (
int i = 0; i < kStateCount; i++) {
178 REALTYPE scaledEigenValue = gEigenValues[eigenIndex][i] * ((REALTYPE)categoryRates[l]);
179 matrixTmp[i] = exp(scaledEigenValue * ((REALTYPE)edgeLengths[u]));
180 firstDerivTmp[i] = scaledEigenValue * matrixTmp[i];
184 for (
int i = 0; i < kStateCount; i++) {
185 for (
int j = 0; j < kStateCount; j++) {
187 REALTYPE sumD1 = 0.0;
188 for (
int k = 0; k < kStateCount; k++) {
189 sum += gCMatrices[eigenIndex][m] * matrixTmp[k];
190 sumD1 += gCMatrices[eigenIndex][m] * firstDerivTmp[k];
194 transitionMat[n] = sum;
196 transitionMat[n] = 0;
197 firstDerivMat[n] = sumD1;
201 transitionMat[n] = 1.0;
202 firstDerivMat[n] = 0.0;
209 for (
int u = 0; u < count; u++) {
210 REALTYPE* transitionMat = transitionMatrices[probabilityIndices[u]];
211 REALTYPE* firstDerivMat = transitionMatrices[firstDerivativeIndices[u]];
212 REALTYPE* secondDerivMat = transitionMatrices[secondDerivativeIndices[u]];
214 for (
int l = 0; l < kCategoryCount; l++) {
216 for (
int i = 0; i < kStateCount; i++) {
217 REALTYPE scaledEigenValue = gEigenValues[eigenIndex][i] * ((REALTYPE)categoryRates[l]);
218 matrixTmp[i] = exp(scaledEigenValue * ((REALTYPE)edgeLengths[u]));
219 firstDerivTmp[i] = scaledEigenValue * matrixTmp[i];
220 secondDerivTmp[i] = scaledEigenValue * firstDerivTmp[i];
224 for (
int i = 0; i < kStateCount; i++) {
225 for (
int j = 0; j < kStateCount; j++) {
227 REALTYPE sumD1 = 0.0;
228 REALTYPE sumD2 = 0.0;
229 for (
int k = 0; k < kStateCount; k++) {
230 sum += gCMatrices[eigenIndex][m] * matrixTmp[k];
231 sumD1 += gCMatrices[eigenIndex][m] * firstDerivTmp[k];
232 sumD2 += gCMatrices[eigenIndex][m] * secondDerivTmp[k];
236 transitionMat[n] = sum;
238 transitionMat[n] = 0;
239 firstDerivMat[n] = sumD1;
240 secondDerivMat[n] = sumD2;
244 transitionMat[n] = 1.0;
245 firstDerivMat[n] = 0.0;
246 secondDerivMat[n] = 0.0;
258 #endif // _EigenDecompositionCube_hpp_