1 /* 2 3 PhyML: a program that computes maximum likelihood phylogenies from 4 DNA or AA homologous sequences. 5 6 Copyright (C) Stephane Guindon. Oct 2003 onward. 7 8 All parts of the source except where indicated are distributed under 9 the GNU public licence. See http://www.opensource.org for details. 10 11 */ 12 13 #include <config.h> 14 15 #ifndef RATES_H 16 #define RATES_H 17 18 #include "utilities.h" 19 #include "spr.h" 20 #include "lk.h" 21 #include "optimiz.h" 22 #include "bionj.h" 23 #include "models.h" 24 #include "free.h" 25 #include "help.h" 26 #include "simu.h" 27 #include "eigen.h" 28 #include "pars.h" 29 #include "alrt.h" 30 #include "time.h" 31 #include "m4.h" 32 #include "draw.h" 33 #include "mcmc.h" 34 #include "stats.h" 35 36 void RATES_Monte_Carlo_Mean_Rates(t_tree *tree); 37 void RATES_Monte_Carlo_Mean_Rates_Pre(t_node *a, t_node *d, t_edge *b, phydbl curr_rate, t_tree *tree); 38 void RATES_Print_Rates(t_tree *tree); 39 void RATES_Print_Rates_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree); 40 t_rate *RATES_Make_Rate_Struct(int n_otu); 41 void RATES_Init_Rate_Struct(t_rate *rates, t_rate *existing_rates, int n_otu); 42 void RATES_Classify_Branches(t_tree *tree); 43 void RATES_Adjust_Rates(t_tree *tree); 44 void RATES_Adjust_Rates_Local_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree); 45 void RATES_Adjust_Rates_Local(t_node *a, t_node *d, t_edge *b1, t_tree *tree); 46 void RATES_Record_T(t_tree *tree); 47 void RATES_Restore_T(t_tree *tree); 48 void RATES_Monte_Carlo_Mean_Rates_Core(phydbl t_lim_sup, phydbl t_lim_inf, phydbl *curr_rate, phydbl *mean_rate, phydbl lexp, phydbl alpha); 49 phydbl RATES_Lk_Rates(t_tree *tree); 50 void RATES_Lk_Rates_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree); 51 void RATES_Fill_Node_Rates_Pre(t_node *a, t_node *d, t_edge *b, phydbl *node_r, t_tree *tree); 52 void RATES_Fill_Node_Rates(phydbl *node_r, t_tree *tree); 53 void RATES_Optimize_Node_Times_Serie_Fixed_Br_Len(t_node *a, t_node *d, t_tree *tree); 54 void RATES_Optimize_Lexp(t_tree *tree); 55 void RATES_Round_Optimize(t_tree *tree); 56 void RATES_Optimize_Lexp(t_tree *tree); 57 void RATES_Optimize_Alpha(t_tree *tree); 58 phydbl RATES_Dmu(phydbl mu, int n_jumps, phydbl dt, phydbl a, phydbl b, phydbl lexp, int min_n, int jps_dens); 59 phydbl RATES_Dr_X_Dx(phydbl r, phydbl mu, phydbl y, phydbl dt, phydbl a, phydbl b, phydbl lexp); 60 phydbl RATES_Dmu_Given_Y_Trpzd(phydbl mu, phydbl y, phydbl dt, phydbl a, phydbl b, phydbl lexp, 61 int nsteps, phydbl beg, phydbl end, phydbl prevs); 62 phydbl RATES_Dmu_Given_Y_Std(phydbl mu, phydbl y, phydbl dt, phydbl a, phydbl b, phydbl lexp); 63 phydbl RATES_Dmu_Given_Y_Romb(phydbl mu, phydbl y, phydbl dt, phydbl a, phydbl b, phydbl lexp); 64 phydbl RATES_Dmu_Given_Y(phydbl mu, phydbl y, phydbl dt, phydbl a, phydbl b, phydbl lexp); 65 phydbl RATES_Dy_Given_Mu(phydbl mu, phydbl y, phydbl dt, phydbl a, phydbl b, phydbl lexp); 66 phydbl RATES_Dmu2_Given_Y_X_Dy_Given_Mu1(phydbl mu1, phydbl mu2, phydbl y, phydbl dt1, phydbl dt2, phydbl a, phydbl b, phydbl lexp); 67 phydbl RATES_Dmu2_Given_Mu1_Trpz(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, phydbl a, phydbl b, phydbl lexp, 68 int nsteps, phydbl beg, phydbl end, phydbl prevs); 69 phydbl RATES_Dmu2_And_Mu1(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, phydbl a, phydbl b, phydbl lexp); 70 phydbl RATES_Dmu2_Given_Mu1(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, phydbl a, phydbl b, phydbl lexp); 71 phydbl RATES_Dmu2_Given_Mu1_Romb(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, phydbl a, phydbl b, phydbl lexp); 72 void RATES_Random_Branch_Lengths(t_tree *tree); 73 void RATES_Bracket_N_Jumps(int *up, int *down, phydbl param); 74 void RATES_Set_Node_Times(t_tree *tree); 75 void RATES_Set_Node_Times_Pre(t_node *a, t_node *d, t_tree *tree); 76 void RATES_Optimize_Node_Times(t_tree *tree); 77 phydbl RATES_Exp_Y(phydbl mu1, phydbl mu2, phydbl dt1, phydbl lexp); 78 phydbl RATES_Dmu2_Given_Mu1_Bis(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, phydbl alpha, phydbl beta, phydbl lexp); 79 void RATES_Replace_Br_Lengths_By_Rates_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree); 80 void RATES_Replace_Br_Lengths_By_Rates(t_tree *tree); 81 82 void RATES_Get_Mean_Rates(t_tree *tree); 83 void RATES_Get_Mean_Rates_Pre(t_node *a, t_node *d, t_edge *b, phydbl r_a, t_tree *tree); 84 void RATES_Expect_Number_Subst(phydbl t_beg, phydbl t_end, phydbl r_beg, int *n_jumps, phydbl *mean_r, phydbl *r_end, t_rate *rates, t_tree *tree); 85 void RATES_Optimize_Clock_Rate(t_tree *tree); 86 phydbl RATES_Dmu1_Given_Lbda_And_Mu2(phydbl lbda, phydbl mu1, phydbl mu2, phydbl alpha, phydbl beta); 87 phydbl RATES_Dmu1_And_Mu2_One_Jump_Trpz(phydbl mu1, phydbl mu2, phydbl a, phydbl b, 88 int nsteps, phydbl beg, phydbl end, phydbl prevs); 89 phydbl RATES_Dmu1_And_Mu2_One_Jump_One_Interval(phydbl mu1, phydbl mu2, phydbl a, phydbl b); 90 phydbl RATES_Dmu1_And_Mu2_One_Jump_Two_Intervals(phydbl dt1, phydbl dt2, phydbl mu1, phydbl mu2, phydbl a, phydbl b); 91 92 phydbl RATES_Dmu1_And_Mu2_One_Jump_Old(phydbl mu1, phydbl mu2, phydbl a, phydbl b); 93 94 phydbl RATES_Dmu2_And_Min_N_Given_Mu1(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, int n_min, phydbl a, phydbl b, phydbl lexp); 95 phydbl RATES_Dmu2_And_Mu1_Given_N(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, int n, phydbl a, phydbl b, phydbl lexp); 96 phydbl RATES_Lk_Rates_Core(phydbl br_r_a, phydbl br_r_d, phydbl nd_r_a, phydbl nd_r_d, int n_a, int n_d, phydbl dt_a, phydbl dt_d, t_tree *tree); 97 void RATES_Init_Triplets(t_tree *tree); 98 phydbl RATES_Lk_Change_One_Time(t_node *n, phydbl new_t, t_tree *tree); 99 void RATES_Update_Triplet(t_node *n, t_tree *tree); 100 void RATES_Print_Triplets(t_tree *tree); 101 phydbl RATES_Lk_Change_One_Rate(t_node *d, phydbl new_rate, t_tree *tree); 102 phydbl RATES_Dmu2_And_Mu1_Given_Min_N(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, int n_min, phydbl a, phydbl b, phydbl lexp); 103 phydbl RATES_Dmu2_And_Mu1_Given_N_Normal(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, int n, phydbl a, phydbl b, phydbl lexp); 104 phydbl RATES_Coeff_Corr(phydbl alpha, phydbl beta, int n1, int n2); 105 phydbl RATES_Dmu2_And_Mu1_Given_N_Full(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, int n, phydbl a, phydbl b, phydbl lexp); 106 phydbl RATES_Dmu1_Given_V_And_N(phydbl mu1, phydbl v, int n, phydbl dt1, phydbl a, phydbl b); 107 phydbl RATES_Yule(t_tree *tree); 108 phydbl RATES_Check_Mean_Rates(t_tree *tree); 109 void RATES_Check_Mean_Rates_Pre(t_node *a, t_node *d, t_edge *b, phydbl *sum, t_tree *tree); 110 void RATES_Discretize_Rates(t_tree *tree); 111 void RATES_Discretize_Rates_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree); 112 phydbl RATES_Dmu_Given_V_And_MinN(phydbl mu, phydbl dt, phydbl v, int minn, phydbl a, phydbl b, phydbl lexp); 113 phydbl RATES_Dmu_One(phydbl mu, phydbl dt, phydbl a, phydbl b, phydbl lexp); 114 phydbl RATES_Compound_Core(phydbl mu1, phydbl mu2, int n1, int n2, phydbl dt1, phydbl dt2, phydbl alpha, phydbl beta, phydbl lexp, phydbl eps, int approx); 115 void RATES_Record_Rates(t_tree *tree); 116 void RATES_Reset_Rates(t_tree *tree); 117 void RATES_Record_Times(t_tree *tree); 118 void RATES_Reset_Times(t_tree *tree); 119 void RATES_Update_T_Rates_Pre(t_node *a, t_node *d, t_tree *tree); 120 void RATES_Update_T_Rates(t_tree *tree); 121 void RATES_Get_Rates_From_Bl(t_tree *tree); 122 phydbl RATES_Compound_Core_Joint(phydbl mu1, phydbl mu2, int n1, int n2, phydbl dt1, phydbl dt2, 123 phydbl alpha, phydbl beta, phydbl lexp, phydbl eps, int approx); 124 phydbl RATES_Dmu_Joint(phydbl mu, int n, phydbl dt, phydbl a, phydbl b, phydbl lexp, int min_n); 125 phydbl RATES_Compound_Core_Marginal(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, phydbl alpha, 126 phydbl beta, phydbl lexp, phydbl eps, int approx); 127 phydbl RATES_Lk_Jumps(t_tree *tree); 128 void RATES_Posterior_Rates(t_tree *tree); 129 void RATES_Posterior_One_Rate(t_node *a, t_node *d, int traversal, t_tree *tree); 130 void RATES_Free_Rates(t_rate *rates); 131 void RATES_Initialize_True_Rates(t_tree *tree); 132 void RATES_Posterior_Times(t_tree *tree); 133 void RATES_Posterior_One_Time(t_node *a, t_node *d, int traversal, t_tree *tree); 134 void RATES_Update_Cur_Bl(t_tree *tree); 135 void RATES_Update_Cur_Bl_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree); 136 void RATES_Get_Cov_Matrix_Rooted(phydbl *unroot_cov, t_tree *tree); 137 void RATES_Get_Cov_Matrix_Rooted_Pre(t_node *a, t_node *d, t_edge *b, phydbl *cov, t_tree *tree); 138 void RATES_Bl_To_Ml(t_tree *tree); 139 void RATES_Bl_To_Ml_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree); 140 int RATES_Check_Node_Times(t_tree *tree); 141 void RATES_Check_Node_Times_Pre(t_node *a, t_node *d, int *err, t_tree *tree); 142 void RATES_Covariance_Mu(t_tree *tree); 143 void RATES_Variance_Mu_Pre(t_node *a, t_node *d, t_tree *tree); 144 void RATES_Fill_Lca_Table(t_tree *tree); 145 void RATES_Posterior_Clock_Rate(t_tree *tree); 146 void RATES_Get_Conditional_Variances(t_tree *tree); 147 void RATES_Get_All_Reg_Coeff(t_tree *tree); 148 void RATES_Posterior_Time_Root(t_tree *tree); 149 void RATES_Get_Trip_Conditional_Variances(t_tree *tree); 150 void RATES_Get_All_Trip_Reg_Coeff(t_tree *tree); 151 void RATES_Check_Lk_Rates(t_tree *tree, int *err); 152 phydbl RATES_Expected_Tree_Length(t_tree *tree); 153 void RATES_Expected_Tree_Length_Pre(t_node *a, t_node *d, phydbl eranc, phydbl *mean, int *n, t_tree *tree); 154 void RATES_Normalise_Rates(t_tree *tree); 155 phydbl RATES_Check_Mean_Rates_True(t_tree *tree); 156 phydbl RATES_Find_Min_Dt_Bisec(phydbl r0, phydbl r1, phydbl t0, phydbl t1, phydbl nu, phydbl threshp, int inf); 157 phydbl RATES_Find_Max_Dt_Bisec(phydbl r0, phydbl r1, phydbl t0, phydbl t1, phydbl nu, phydbl threshp, int inf); 158 void RATES_Min_Max_Interval(phydbl u0, phydbl u1, phydbl u2, phydbl u3, phydbl t0, phydbl t2, phydbl t3, 159 phydbl *t_min, phydbl *t_max, phydbl nu, phydbl p_thresh, t_tree *tree); 160 161 162 phydbl RATES_Get_Correction_Factor(phydbl mode, phydbl sd, int *err, t_tree *tree); 163 phydbl RATES_Average_Substitution_Rate(t_tree *tree); 164 void RATES_Update_Norm_Fact(t_tree *tree); 165 void RATES_Update_Mean_Br_Len(int iter, t_tree *tree); 166 void RATES_Update_Cov_Br_Len(int iter, t_tree *tree); 167 void RATES_Set_Mean_L(t_tree *tree); 168 void RATES_Fill_All_Param(t_rate *rate, t_tree *tree); 169 void RATES_Record_Rates(t_tree *tree); 170 void RATES_Reset_Rates(t_tree *tree); 171 phydbl RATES_Average_Rate(t_tree *tree); 172 void RATES_Set_Clock_And_Nu_Max(t_tree *tree); 173 phydbl RATES_Lk_Linreg(t_tree *tree); 174 phydbl RATES_Get_Mean_Rate_In_Subtree(t_node *root, t_tree *tree); 175 void RATES_Get_Mean_Rate_In_Subtree_Pre(t_node *a, t_node *d, phydbl *sum, int *n, t_tree *tree); 176 char *RATES_Get_Model_Name(int model); 177 void RATES_Get_Survival_Ranks(t_tree *tree); 178 void RATES_Bl_To_Bl(t_tree *tree); 179 void RATES_Bl_To_Bl_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree); 180 void RATES_Set_Birth_Rate_Boundaries(t_tree *tree); 181 void RATES_Copy_Rate_Struct(t_rate *from, t_rate *to, int n_otu); 182 void RATES_Duplicate_Calib_Struct(t_tree *from, t_tree *to); 183 int RATES_Check_Edge_Length_Consistency(t_tree *tree); 184 phydbl RATES_Realized_Substitution_Rate(t_tree *tree); 185 186 #endif 187