7 #ifndef _EigenDecompositionSquare_hpp_
8 #define _EigenDecompositionSquare_hpp_
9 #include "EigenDecompositionSquare.h"
10 #include "libhmsbeagle/beagle.h"
23 BEAGLE_CPU_EIGEN_TEMPLATE
24 EigenDecompositionSquare<BEAGLE_CPU_EIGEN_GENERIC>::EigenDecompositionSquare(
int decompositionCount,
28 : EigenDecomposition<BEAGLE_CPU_EIGEN_GENERIC>(decompositionCount,stateCount,categoryCount, flags) {
30 isComplex = kFlags & BEAGLE_FLAG_EIGEN_COMPLEX;
33 kEigenValuesSize = 2 * kStateCount;
35 kEigenValuesSize = kStateCount;
37 this->gEigenValues = (REALTYPE**) malloc(
sizeof(REALTYPE*) * kEigenDecompCount);
38 if (gEigenValues == NULL)
39 throw std::bad_alloc();
41 gEMatrices = (REALTYPE**) malloc(
sizeof(REALTYPE*) * kEigenDecompCount);
42 if (gEMatrices == NULL)
43 throw std::bad_alloc();
45 gIMatrices = (REALTYPE**) malloc(
sizeof(REALTYPE*) * kEigenDecompCount);
46 if (gIMatrices == NULL)
47 throw std::bad_alloc();
49 for (
int i = 0; i < kEigenDecompCount; i++) {
50 gEMatrices[i] = (REALTYPE*) malloc(
sizeof(REALTYPE) * kStateCount * kStateCount);
51 if (gEMatrices[i] == NULL)
52 throw std::bad_alloc();
54 gIMatrices[i] = (REALTYPE*) malloc(
sizeof(REALTYPE) * kStateCount * kStateCount);
55 if (gIMatrices[i] == NULL)
56 throw std::bad_alloc();
58 gEigenValues[i] = (REALTYPE*) malloc(
sizeof(REALTYPE) * kEigenValuesSize);
59 if (gEigenValues[i] == NULL)
60 throw std::bad_alloc();
63 matrixTmp = (REALTYPE*) malloc(
sizeof(REALTYPE) * kStateCount * kStateCount);
66 BEAGLE_CPU_EIGEN_TEMPLATE
67 EigenDecompositionSquare<BEAGLE_CPU_EIGEN_GENERIC>::~EigenDecompositionSquare() {
69 for(
int i=0; i<kEigenDecompCount; i++) {
72 free(gEigenValues[i]);
83 template<
typename REALTYPE>
84 void transposeSquareMatrix(REALTYPE* mat,
86 for (
int i = 0; i < size - 1; i++) {
87 for (
int j = i + 1; j < size; j++) {
88 REALTYPE tmp = mat[i * size + j];
89 mat[i * size + j] = mat[j * size + i];
90 mat[j * size + i] = tmp;
95 BEAGLE_CPU_EIGEN_TEMPLATE
96 void EigenDecompositionSquare<BEAGLE_CPU_EIGEN_GENERIC>::setEigenDecomposition(
int eigenIndex,
97 const double* inEigenVectors,
98 const double* inInverseEigenVectors,
99 const double* inEigenValues) {
101 beagleMemCpy(gEigenValues[eigenIndex],inEigenValues,kEigenValuesSize);
102 const int len = kStateCount * kStateCount;
103 beagleMemCpy(gEMatrices[eigenIndex],inEigenVectors,len);
104 beagleMemCpy(gIMatrices[eigenIndex],inInverseEigenVectors,len);
105 if (kFlags & BEAGLE_FLAG_INVEVEC_TRANSPOSED)
106 transposeSquareMatrix(gIMatrices[eigenIndex], kStateCount);
109 BEAGLE_CPU_EIGEN_TEMPLATE
110 void EigenDecompositionSquare<BEAGLE_CPU_EIGEN_GENERIC>::updateTransitionMatrices(
int eigenIndex,
111 const int* probabilityIndices,
112 const int* firstDerivativeIndices,
113 const int* secondDerivativeIndices,
114 const double* edgeLengths,
115 const double* categoryRates,
116 REALTYPE** transitionMatrices,
119 const REALTYPE* Ievc = gIMatrices[eigenIndex];
120 const REALTYPE* Evec = gEMatrices[eigenIndex];
121 const REALTYPE* Eval = gEigenValues[eigenIndex];
122 const REALTYPE* EvalImag = Eval + kStateCount;
123 for (
int u = 0; u < count; u++) {
124 REALTYPE* transitionMat = transitionMatrices[probabilityIndices[u]];
125 const double edgeLength = edgeLengths[u];
127 for (
int l = 0; l < kCategoryCount; l++) {
128 const REALTYPE distance = categoryRates[l] * edgeLength;
129 for(
int i=0; i<kStateCount; i++) {
130 if (!isComplex || EvalImag[i] == 0) {
131 const REALTYPE tmp = exp(Eval[i] * distance);
132 for(
int j=0; j<kStateCount; j++) {
133 matrixTmp[i*kStateCount+j] = Ievc[i*kStateCount+j] * tmp;
138 const REALTYPE b = EvalImag[i];
139 const REALTYPE expat = exp(Eval[i] * distance);
140 const REALTYPE expatcosbt = expat * cos(b * distance);
141 const REALTYPE expatsinbt = expat * sin(b * distance);
142 for(
int j=0; j<kStateCount; j++) {
143 matrixTmp[ i*kStateCount+j] = expatcosbt * Ievc[ i*kStateCount+j] +
144 expatsinbt * Ievc[i2*kStateCount+j];
145 matrixTmp[i2*kStateCount+j] = expatcosbt * Ievc[i2*kStateCount+j] -
146 expatsinbt * Ievc[ i*kStateCount+j];
154 for(
int i=0; i<16; i++)
155 fprintf(stderr,
" %7.5e,",matrixTmp[i]);
156 fprintf(stderr,
"] -- complex debug\n");
161 for (
int i = 0; i < kStateCount; i++) {
162 for (
int j = 0; j < kStateCount; j++) {
164 for (
int k = 0; k < kStateCount; k++)
165 sum += Evec[i*kStateCount+k] * matrixTmp[k*kStateCount+j];
167 transitionMat[n] = sum;
169 transitionMat[n] = 0;
173 transitionMat[n] = 1.0;
179 if (DEBUGGING_OUTPUT) {
180 int kMatrixSize = kStateCount * kStateCount;
181 fprintf(stderr,
"transitionMat index=%d brlen=%.5f\n", probabilityIndices[u], edgeLengths[u]);
182 for (
int w = 0; w < (20 > kMatrixSize ? 20 : kMatrixSize); ++w)
183 fprintf(stderr,
"transitionMat[%d] = %.5f\n", w, transitionMat[w]);