HMSBEAGLE  1.0.0
EigenDecompositionCube.hpp
1 /*
2  * EigenDecompositionCube.cpp
3  *
4  * Created on: Sep 24, 2009
5  * Author: msuchard
6  */
7 #ifndef _EigenDecompositionCube_hpp_
8 #define _EigenDecompositionCube_hpp_
9 
10 #include "libhmsbeagle/CPU/EigenDecompositionCube.h"
11 
12 
13 namespace beagle {
14 namespace cpu {
15 
16 #if defined (BEAGLE_IMPL_DEBUGGING_OUTPUT) && BEAGLE_IMPL_DEBUGGING_OUTPUT
17 const bool DEBUGGING_OUTPUT = true;
18 #else
19 const bool DEBUGGING_OUTPUT = false;
20 #endif
21 
22 BEAGLE_CPU_EIGEN_TEMPLATE
23 EigenDecompositionCube<BEAGLE_CPU_EIGEN_GENERIC>::EigenDecompositionCube(int decompositionCount,
24  int stateCount,
25  int categoryCount,
26  long flags)
27  : EigenDecomposition<BEAGLE_CPU_EIGEN_GENERIC>(decompositionCount,
28  stateCount,
29  categoryCount,
30  flags) {
31  gEigenValues = (REALTYPE**) malloc(sizeof(REALTYPE*) * kEigenDecompCount);
32  if (gEigenValues == NULL)
33  throw std::bad_alloc();
34 
35  gCMatrices = (REALTYPE**) malloc(sizeof(REALTYPE*) * kEigenDecompCount);
36  if (gCMatrices == NULL)
37  throw std::bad_alloc();
38 
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();
43 
44  gEigenValues[i] = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount);
45  if (gEigenValues[i] == NULL)
46  throw std::bad_alloc();
47  }
48 
49  matrixTmp = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount);
50  firstDerivTmp = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount);
51  secondDerivTmp = (REALTYPE*) malloc(sizeof(REALTYPE) * kStateCount);
52 }
53 
54 BEAGLE_CPU_EIGEN_TEMPLATE
55 EigenDecompositionCube<BEAGLE_CPU_EIGEN_GENERIC>::~EigenDecompositionCube() {
56 
57  for(int i=0; i<kEigenDecompCount; i++) {
58  free(gCMatrices[i]);
59  free(gEigenValues[i]);
60  }
61  free(gCMatrices);
62  free(gEigenValues);
63  free(matrixTmp);
64  free(firstDerivTmp);
65  free(secondDerivTmp);
66 }
67 
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) {
73 
74  if (kFlags & BEAGLE_FLAG_INVEVEC_STANDARD) {
75  int l = 0;
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];
82  l++;
83  }
84  }
85  }
86  } else {
87  int l = 0;
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)];
94  l++;
95  }
96  }
97  }
98  }
99 
100 }
101 
102 #define UNROLL
103 
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,
112  int count) {
113 #ifdef UNROLL
114  int stateCountModFour = (kStateCount / 4) * 4;
115 #endif
116 
117  if (firstDerivativeIndices == NULL && secondDerivativeIndices == NULL) {
118  for (int u = 0; u < count; u++) {
119  REALTYPE* transitionMat = transitionMatrices[probabilityIndices[u]];
120  int n = 0;
121  for (int l = 0; l < kCategoryCount; l++) {
122 
123  for (int i = 0; i < kStateCount; i++) {
124  matrixTmp[i] = exp(gEigenValues[eigenIndex][i] * ((REALTYPE)edgeLengths[u] * categoryRates[l]));
125  }
126 
127  REALTYPE* tmpCMatrices = gCMatrices[eigenIndex];
128  for (int i = 0; i < kStateCount; i++) {
129  for (int j = 0; j < kStateCount; j++) {
130  REALTYPE sum = 0.0;
131 #ifdef UNROLL
132  int k = 0;
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];
138  }
139  for (; k < kStateCount; k++) {
140  sum += tmpCMatrices[k] * matrixTmp[k];
141  }
142  tmpCMatrices += kStateCount;
143 #else
144  for (int k = 0; k < kStateCount; k++) {
145  sum += *tmpCMatrices++ * matrixTmp[k];
146  }
147 #endif
148  if (sum > 0)
149  transitionMat[n] = sum;
150  else
151  transitionMat[n] = 0;
152  n++;
153  }
154 if (T_PAD != 0) {
155  transitionMat[n] = 1.0;
156  n += T_PAD;
157 }
158  }
159  }
160 
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]);
166  }
167  }
168 
169 
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]];
174  int n = 0;
175  for (int l = 0; l < kCategoryCount; l++) {
176 
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];
181  }
182 
183  int m = 0;
184  for (int i = 0; i < kStateCount; i++) {
185  for (int j = 0; j < kStateCount; j++) {
186  REALTYPE sum = 0.0;
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];
191  m++;
192  }
193  if (sum > 0)
194  transitionMat[n] = sum;
195  else
196  transitionMat[n] = 0;
197  firstDerivMat[n] = sumD1;
198  n++;
199  }
200 if (T_PAD != 0) {
201  transitionMat[n] = 1.0;
202  firstDerivMat[n] = 0.0;
203  n += T_PAD;
204 }
205  }
206  }
207  }
208  } else {
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]];
213  int n = 0;
214  for (int l = 0; l < kCategoryCount; l++) {
215 
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];
221  }
222 
223  int m = 0;
224  for (int i = 0; i < kStateCount; i++) {
225  for (int j = 0; j < kStateCount; j++) {
226  REALTYPE sum = 0.0;
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];
233  m++;
234  }
235  if (sum > 0)
236  transitionMat[n] = sum;
237  else
238  transitionMat[n] = 0;
239  firstDerivMat[n] = sumD1;
240  secondDerivMat[n] = sumD2;
241  n++;
242  }
243 if (T_PAD != 0) {
244  transitionMat[n] = 1.0;
245  firstDerivMat[n] = 0.0;
246  secondDerivMat[n] = 0.0;
247  n += T_PAD;
248 }
249  }
250  }
251  }
252  }
253 }
254 
255 } // cpu
256 } // beagle
257 
258 #endif // _EigenDecompositionCube_hpp_
259