HMSBEAGLE  1.0.0
EigenDecompositionSquare.hpp
1 /*
2  * EigenDecompositionSquare.cpp
3  *
4  * Created on: Sep 24, 2009
5  * Author: msuchard
6  */
7 #ifndef _EigenDecompositionSquare_hpp_
8 #define _EigenDecompositionSquare_hpp_
9 #include "EigenDecompositionSquare.h"
10 #include "libhmsbeagle/beagle.h"
11 
12 //#if defined (BEAGLE_IMPL_DEBUGGING_OUTPUT) && BEAGLE_IMPL_DEBUGGING_OUTPUT
13 //const bool DEBUGGING_OUTPUT = true;
14 //#else
15 //const bool DEBUGGING_OUTPUT = false;
16 //#endif
17 
18 //#define DEBUG_COMPLEX
19 
20 namespace beagle {
21 namespace cpu {
22 
23 BEAGLE_CPU_EIGEN_TEMPLATE
24 EigenDecompositionSquare<BEAGLE_CPU_EIGEN_GENERIC>::EigenDecompositionSquare(int decompositionCount,
25  int stateCount,
26  int categoryCount,
27  long flags)
28  : EigenDecomposition<BEAGLE_CPU_EIGEN_GENERIC>(decompositionCount,stateCount,categoryCount, flags) {
29 
30  isComplex = kFlags & BEAGLE_FLAG_EIGEN_COMPLEX;
31 
32  if (isComplex)
33  kEigenValuesSize = 2 * kStateCount;
34  else
35  kEigenValuesSize = kStateCount;
36 
37  this->gEigenValues = (REALTYPE**) malloc(sizeof(REALTYPE*) * kEigenDecompCount);
38  if (gEigenValues == NULL)
39  throw std::bad_alloc();
40 
41  gEMatrices = (REALTYPE**) malloc(sizeof(REALTYPE*) * kEigenDecompCount);
42  if (gEMatrices == NULL)
43  throw std::bad_alloc();
44 
45  gIMatrices = (REALTYPE**) malloc(sizeof(REALTYPE*) * kEigenDecompCount);
46  if (gIMatrices == NULL)
47  throw std::bad_alloc();
48 
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();
53 
54  gIMatrices[i] = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount * kStateCount);
55  if (gIMatrices[i] == NULL)
56  throw std::bad_alloc();
57 
58  gEigenValues[i] = (REALTYPE*) malloc(sizeof(REALTYPE) * kEigenValuesSize);
59  if (gEigenValues[i] == NULL)
60  throw std::bad_alloc();
61  }
62 
63  matrixTmp = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount * kStateCount);
64 }
65 
66 BEAGLE_CPU_EIGEN_TEMPLATE
67 EigenDecompositionSquare<BEAGLE_CPU_EIGEN_GENERIC>::~EigenDecompositionSquare() {
68 
69  for(int i=0; i<kEigenDecompCount; i++) {
70  free(gEMatrices[i]);
71  free(gIMatrices[i]);
72  free(gEigenValues[i]);
73  }
74  free(gEMatrices);
75  free(gIMatrices);
76  free(gEigenValues);
77  free(matrixTmp);
78 }
79 
83 template<typename REALTYPE>
84 void transposeSquareMatrix(REALTYPE* mat,
85  int size) {
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;
91  }
92  }
93 }
94 
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) {
100 
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) // TODO: optimize, might not need to transpose here
106  transposeSquareMatrix(gIMatrices[eigenIndex], kStateCount);
107 }
108 
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,
117  int count) {
118 
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];
126  int n = 0;
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;
134  }
135  } else {
136  // 2 x 2 conjugate block
137  int i2 = i + 1;
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];
147  }
148  i++; // processed two conjugate rows
149  }
150  }
151 
152 #ifdef DEBUG_COMPLEX
153  fprintf(stderr,"[");
154  for(int i=0; i<16; i++)
155  fprintf(stderr," %7.5e,",matrixTmp[i]);
156  fprintf(stderr,"] -- complex debug\n");
157  exit(0);
158 #endif
159 
160 
161  for (int i = 0; i < kStateCount; i++) {
162  for (int j = 0; j < kStateCount; j++) {
163  REALTYPE sum = 0.0;
164  for (int k = 0; k < kStateCount; k++)
165  sum += Evec[i*kStateCount+k] * matrixTmp[k*kStateCount+j];
166  if (sum > 0)
167  transitionMat[n] = sum;
168  else
169  transitionMat[n] = 0;
170  n++;
171  }
172 if (T_PAD != 0) {
173  transitionMat[n] = 1.0;
174  n += T_PAD;
175 }
176  }
177  }
178 
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]);
184  }
185  }
186 }
187 
188 }
189 }
190 
191 #endif
192