HMSBEAGLE  1.0.0
libhmsbeagle/CPU/EigenDecompositionCube.hpp
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