HMSBEAGLE  1.0.0
libhmsbeagle/CPU/EigenDecomposition.h
00001 /*
00002  * EigenDecomposition.h
00003  *
00004  *  Created on: Sep 24, 2009
00005  *      Author: msuchard
00006  */
00007 
00008 #ifndef EIGENDECOMPOSITION_H_
00009 #define EIGENDECOMPOSITION_H_
00010 
00011 #include <cstdio>
00012 #include <cstdlib>
00013 #include <iostream>
00014 #include <cstring>
00015 #include <cmath>
00016 #include <cassert>
00017 #include <vector>
00018 
00019 #define BEAGLE_CPU_EIGEN_GENERIC        REALTYPE, T_PAD
00020 #define BEAGLE_CPU_EIGEN_TEMPLATE       template <typename REALTYPE, int T_PAD>
00021 
00022 namespace beagle {
00023 namespace cpu {
00024 
00025 BEAGLE_CPU_EIGEN_TEMPLATE
00026 class EigenDecomposition {
00027         
00028 protected:
00029     REALTYPE** gEigenValues;
00030     int kStateCount;
00031     int kEigenDecompCount;
00032     int kCategoryCount;
00033         long kFlags;
00034     REALTYPE* matrixTmp;
00035     REALTYPE* firstDerivTmp;
00036     REALTYPE* secondDerivTmp;
00037     
00038 public:
00039         EigenDecomposition(int decompositionCount,
00040                                            int stateCount,
00041                                            int categoryCount,
00042                        long flags)
00043                                            {
00044 
00045                                                         kEigenDecompCount = decompositionCount;
00046                                                         kStateCount = stateCount;
00047                                                         kCategoryCount = categoryCount;
00048                             kFlags = flags;
00049                                                 };
00050         
00051         virtual ~EigenDecomposition() {};
00052         
00053     // sets the Eigen decomposition for a given matrix
00054     //
00055     // matrixIndex the matrix index to update
00056     // eigenVectors an array containing the Eigen Vectors
00057     // inverseEigenVectors an array containing the inverse Eigen Vectors
00058     // eigenValues an array containing the Eigen Values
00059     virtual void setEigenDecomposition(int eigenIndex,
00060                               const double* inEigenVectors,
00061                               const double* inInverseEigenVectors,
00062                               const double* inEigenValues) = 0;
00063                 
00064     // calculate a transition probability matrices for a given list of node. This will
00065     // calculate for all categories (and all matrices if more than one is being used).
00066     //
00067     // nodeIndices an array of node indices that require transition probability matrices
00068     // edgeLengths an array of expected lengths in substitutions per site
00069     // count the number of elements in the above arrays
00070     virtual void updateTransitionMatrices(int eigenIndex,
00071                                  const int* probabilityIndices,
00072                                  const int* firstDerivativeIndices,
00073                                  const int* secondDerivativeIndices,
00074                                  const double* edgeLengths,
00075                                  const double* categoryRates,
00076                                  REALTYPE** transitionMatrices,
00077                                  int count) = 0;
00078 
00079 };
00080 
00081 }
00082 }
00083 
00084 #endif /* EIGENDECOMPOSITION_H_ */