1 /* 2 * author: Imran Fanaswala 3 */ 4 5 #ifndef BEAGLE_UTILS_H 6 #define BEAGLE_UTILS_H 7 8 #include "assert.h" 9 #include "libhmsbeagle/beagle.h" 10 #include "utilities.h" 11 #define UNINITIALIZED -42 12 13 /* 14 BEAGLE is a library that abstracts the core computations common to most 15 phylogenetic packages. The following paragraphs serve as a "design document" for 16 the implementation details of how BEAGLE has been integrated with PhyML. The 17 reader is expected to be familiar with phylogenetics, skimmed the BEAGLE API, 18 and seen a couple of its examples. 19 20 PARTIAL LIKELIHOODS: 21 - For each internal edge, PhyML maintains a vector of partials that lie to the 22 "right" and "left", to faciliate both post-order and pre-order tree traversal. 23 Therefore each t_edge has a `p_lk_right` and `p_lk_left` field 24 - PhyML represents a tip as a vector of partials (where one bit is "1", and 25 the rest are "0"). Therefore, a t_edge connected to a tip has `p_lk_tip` field. 26 Also by convention, the tip always lies to the right of the edge its connected 27 to. In other words, if b->right->tax is True, then b->p_lk_tip holds the tip 28 partials. 29 - Therefore for BEAGLE, we allocate N+(2*num_edges) partial buffers. 30 - Recall that BEAGLE requires each partial buffer to be uniquely indexed. Thus 31 each edge also has an additional `p_lk_rght_idx` and `p_lk_left_idx` field. 32 These fields are used to create a BeagleOperation struct during the call to 33 beagleUpdatePartials(). 34 - In PhyML, Update_P_Lk() updates the partials for a specific edge lieing along 35 a specific "direction". The "direction" depends on which tree-traversal is being 36 employed. In other words, the tree-traversal dictates whether we are computing 37 "up"/"right" or "down"/"left" partials relative to the t_node*. Specifically, 38 Set_All_P_Lk() handles the nitty gritty of determining whether we are computing 39 "right" or "left" partials. Therefore this function has been modified to also 40 set the appropriate BEAGLE indices. 41 42 P-MATS and EDGELIKELIHOODS: 43 - (1) Observe that PhyML does eigen decomposition for only GTR and AA models and 44 thus the P-matrix is computed using Beagle. This implies the usage of 45 SetEigenDecomposition() and UpdateTransitionMatricies() calls. 46 - (2) On the other hand, in case of simpler models (JC69, F81, etc), the 47 P-Matrix is explicitly set in BEAGLE thus implying the usage of 48 SetTransitionMatrix(). 49 - The choice between (1) and (2) is made in Update_PMat_At_Given_Edge(). BTW, 50 upon reading Update_PMat_At_Given_Edge(), it may appear that the P-Matrix is 51 being computed in PhyML *and* in BEAGLE; afterall because PMat() function in all 52 cases. However, when BEAGLE is enabled, observe that in PMat() essentially does 53 nothing. 54 - Lk_Core() computes the edgelikelihoods, and calc_edgelks_beagle() is the 55 corresponding BEAGLE function that does the same. 56 57 58 */ 59 60 61 62 #define TODO_BEAGLE "TODO. This codepath has not been implemented in PhyML-X, please post your usecase on the PhyML discussion list" 63 64 int create_beagle_instance(t_tree* tree, int quiet, option* io); 65 int finalize_beagle_instance(t_tree* tree); 66 void update_beagle_partials(t_tree* tree, t_edge* b, t_node* d); 67 void update_beagle_ras(t_mod* mod); 68 void update_beagle_efrqs(t_mod* mod); 69 void update_beagle_eigen(t_mod* mod); 70 void calc_edgelks_beagle(t_edge* b, t_tree* tree); 71 double* int_to_double(const int* src, int num_elems); 72 double* short_to_double(const short* src, int num_elems); 73 double* float_to_double(const phydbl *src, int num_elems); 74 75 76 77 #endif // BEAGLE_UTILS_H 78