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 /* Routines for molecular clock trees and molecular dating */
14
15
16 #include "rates.h"
17
18 #ifdef RWRAPPER
19 #include <R.h>
20 #endif
21
22 //////////////////////////////////////////////////////////////
23 //////////////////////////////////////////////////////////////
24
RATES_Lk_Rates(t_tree * tree)25 phydbl RATES_Lk_Rates(t_tree *tree)
26 {
27
28 if(tree->eval_rlnL == NO) return UNLIKELY;
29
30 tree->rates->c_lnL_rates = .0;
31
32 RATES_Lk_Rates_Pre(tree->n_root,tree->n_root->v[2],NULL,tree);
33 RATES_Lk_Rates_Pre(tree->n_root,tree->n_root->v[1],NULL,tree);
34
35 if(isnan(tree->rates->c_lnL_rates)) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
36
37 return tree->rates->c_lnL_rates;
38 }
39
40 //////////////////////////////////////////////////////////////
41 //////////////////////////////////////////////////////////////
42
RATES_Lk_Rates_Pre(t_node * a,t_node * d,t_edge * b,t_tree * tree)43 void RATES_Lk_Rates_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
44 {
45 int i;
46 phydbl log_dens,mu_a,mu_d,r_a,r_d,dt_a,dt_d;
47 int n_a,n_d;
48
49 log_dens = -1.;
50
51 if(d->anc != a)
52 {
53 PhyML_Fprintf(stderr,"\n. d=%d d->anc=%d a=%d root=%d",d->num,d->anc->num,a->num,tree->n_root->num);
54 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
55 assert(FALSE);
56 }
57
58 dt_a = -1.;
59 if(a != tree->n_root) dt_a = tree->times->nd_t[a->num] - tree->times->nd_t[a->anc->num];
60
61 mu_a = tree->rates->br_r[a->num];
62 r_a = tree->rates->nd_r[a->num];
63 n_a = tree->times->n_jps[a->num];
64
65 dt_d = FABS(tree->times->nd_t[d->num] - tree->times->nd_t[a->num]);
66 mu_d = tree->rates->br_r[d->num];
67 r_d = tree->rates->nd_r[d->num];
68 n_d = tree->times->n_jps[d->num];
69
70 log_dens = RATES_Lk_Rates_Core(mu_a,mu_d,r_a,r_d,n_a,n_d,dt_a,dt_d,tree);
71 tree->rates->c_lnL_rates += log_dens;
72
73 if(isnan(tree->rates->c_lnL_rates))
74 {
75 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
76 MCMC_Print_Param(tree->mcmc,tree);
77 Exit("\n");
78 }
79
80 tree->rates->triplet[a->num] += log_dens;
81
82
83 if(d->tax) return;
84 else
85 {
86 for(i=0;i<3;i++)
87 {
88 if((d->v[i] != a) && (d->b[i] != tree->e_root))
89 {
90 RATES_Lk_Rates_Pre(d,d->v[i],d->b[i],tree);
91 }
92 }
93 }
94 }
95
96 //////////////////////////////////////////////////////////////
97 //////////////////////////////////////////////////////////////
98
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)99 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)
100 {
101 phydbl log_dens;
102 phydbl mean,sd;
103 phydbl min_r, max_r;
104 phydbl cr;
105
106 cr = tree->rates->clock_r;
107 log_dens = UNLIKELY;
108 mean = sd = -1.;
109 min_r = tree->rates->min_rate;
110 max_r = tree->rates->max_rate;
111
112 dt_d = MAX(0.5,dt_d); // We give only one decimal when printing out node heights. It is therefore a fair approximation
113
114 switch(tree->rates->model)
115 {
116 case THORNE :
117 {
118 PhyML_Fprintf(stderr,"\n. Not implemented yet.");
119 assert(FALSE);
120 break;
121 }
122 case GUINDON :
123 {
124 int err;
125
126 min_r = tree->rates->min_rate;
127 max_r = tree->rates->max_rate;
128
129 nd_r_d = log(nd_r_d*cr);
130 nd_r_a = log(nd_r_a*cr);
131 min_r = log(min_r*cr);
132 max_r = log(max_r*cr);
133
134 sd = sqrt(tree->rates->nu*dt_d);
135 mean = nd_r_a - .5*sd*sd;
136
137 // log(density(log(Rate)))
138 log_dens = Log_Dnorm_Trunc(nd_r_d,mean,sd,min_r,max_r,&err);
139 // log(density(Rate))
140 log_dens -= log(exp(nd_r_d)/cr);
141
142 if(err)
143 {
144 PhyML_Fprintf(stderr,"\n. Run: %d",tree->mcmc->run);
145 PhyML_Fprintf(stderr,"\n. br_r_d=%f mean=%f sd=%f min_r=%f max_r=%f dt_d=%f",br_r_d,mean,sd,min_r,max_r,dt_d);
146 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
147 Exit("\n");
148 }
149 break;
150 }
151 case LOGNORMAL :
152 {
153 int err;
154 phydbl log_br_r_d = log(br_r_d);
155 log_dens = Log_Dnorm_Trunc(log_br_r_d,0.0,tree->rates->nu,tree->rates->min_rate,tree->rates->max_rate,&err);
156 log_dens -= log_br_r_d;
157 break;
158 }
159 case STRICTCLOCK :
160 {
161 log_dens = 0.0;
162 break;
163 }
164
165 default :
166 {
167 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
168 Warn_And_Exit("");
169 }
170 }
171
172 if(isnan(log_dens))
173 {
174 PhyML_Fprintf(stderr,"\n. Run=%4d br_r_d=%f br_r_a=%f dt_d=%f dt_a=%f nu=%f log_dens=%G sd=%f mean=%f\n",
175 tree->mcmc->run,
176 br_r_d,br_r_a,dt_d,dt_a,tree->rates->nu,log_dens,
177 sd,mean);
178 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
179 Exit("\n");
180 }
181
182 return log_dens;
183 }
184
185 //////////////////////////////////////////////////////////////
186 //////////////////////////////////////////////////////////////
187
RATES_Lk_Change_One_Rate(t_node * d,phydbl new_rate,t_tree * tree)188 phydbl RATES_Lk_Change_One_Rate(t_node *d, phydbl new_rate, t_tree *tree)
189 {
190 tree->rates->br_r[d->num] = new_rate;
191 RATES_Update_Triplet(d,tree);
192 RATES_Update_Triplet(d->anc,tree);
193 return(tree->rates->c_lnL_rates);
194 }
195
196 //////////////////////////////////////////////////////////////
197 //////////////////////////////////////////////////////////////
198
199
RATES_Lk_Change_One_Time(t_node * n,phydbl new_t,t_tree * tree)200 phydbl RATES_Lk_Change_One_Time(t_node *n, phydbl new_t, t_tree *tree)
201 {
202 if(n == tree->n_root)
203 {
204 PhyML_Fprintf(stderr,"\n. Moving the time of the root t_node is not permitted.");
205 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
206 Warn_And_Exit("");
207 }
208 else
209 {
210 int i;
211
212 tree->times->nd_t[n->num] = new_t;
213
214 RATES_Update_Triplet(n,tree);
215
216 for(i=0;i<3;i++)
217 {
218 if(n->b[i] != tree->e_root) RATES_Update_Triplet(n->v[i],tree);
219 else RATES_Update_Triplet(tree->n_root,tree);
220 }
221 }
222 return(tree->rates->c_lnL_rates);
223 }
224
225 //////////////////////////////////////////////////////////////
226 //////////////////////////////////////////////////////////////
227
RATES_Update_Triplet(t_node * n,t_tree * tree)228 void RATES_Update_Triplet(t_node *n, t_tree *tree)
229 {
230 phydbl curr_triplet,new_triplet;
231 phydbl dt0,dt1,dt2;
232 phydbl mu1_mu0,mu2_mu0;
233 phydbl mu0,mu1,mu2;
234 phydbl r0,r1,r2;
235 int n0,n1,n2;
236 int i;
237 t_node *v1,*v2;
238
239 if(n->tax) return;
240
241 curr_triplet = tree->rates->triplet[n->num];
242
243 dt0 = dt1 = dt2 = -100.0;
244 r0 = r1 = r2 = 0.0;
245
246 if(n == tree->n_root)
247 {
248 phydbl log_dens;
249
250 log_dens = 0.0;
251
252 dt0 = tree->times->nd_t[tree->n_root->v[2]->num] - tree->times->nd_t[tree->n_root->num];
253 dt1 = tree->times->nd_t[tree->n_root->v[1]->num] - tree->times->nd_t[tree->n_root->num];
254
255 mu0 = tree->rates->br_r[tree->n_root->v[2]->num];
256 mu1 = tree->rates->br_r[tree->n_root->v[1]->num];
257
258 r0 = tree->rates->nd_r[tree->n_root->v[2]->num];
259 r1 = tree->rates->nd_r[tree->n_root->v[1]->num];
260
261 n0 = tree->times->n_jps[tree->n_root->v[2]->num];
262 n1 = tree->times->n_jps[tree->n_root->v[1]->num];
263
264 switch(tree->rates->model)
265 {
266 case COMPOUND_COR : case COMPOUND_NOCOR :
267 {
268 log_dens = RATES_Dmu(mu0,n0,dt0,tree->rates->nu,1./tree->rates->nu,tree->rates->lexp,0,1);
269 log_dens *= RATES_Dmu(mu1,n1,dt1,tree->rates->nu,1./tree->rates->nu,tree->rates->lexp,0,1);
270 log_dens = log(log_dens);
271 break;
272 }
273 case EXPONENTIAL :
274 {
275 log_dens = Dexp(mu0,tree->rates->lexp) * Dexp(mu1,tree->rates->lexp);
276 log_dens = log(log_dens);
277 break;
278 }
279 case LOGNORMAL :
280 {
281 log_dens = Dgamma(mu0,tree->rates->nu,1./tree->rates->nu) * Dgamma(mu1,tree->rates->nu,1./tree->rates->nu);
282 log_dens = log(log_dens);
283 break;
284 }
285 case THORNE :
286 {
287 int err;
288 phydbl mean0,sd0;
289 phydbl mean1,sd1;
290
291
292 sd0 = SQRT(tree->rates->nu*dt0);
293 mean0 = 1.0;
294
295 sd1 = SQRT(tree->rates->nu*dt1);
296 mean1 = 1.0;
297
298 log_dens =
299 Log_Dnorm_Trunc(mu0,mean0,sd0,tree->rates->min_rate,tree->rates->max_rate,&err) +
300 Log_Dnorm_Trunc(mu1,mean1,sd1,tree->rates->min_rate,tree->rates->max_rate,&err);
301
302 break;
303 }
304 case GUINDON :
305 {
306 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
307 Exit("\n. Not implemented yet.\n");
308 break;
309 }
310 default :
311 {
312 Exit("\n. Model not implemented yet.\n");
313 break;
314 }
315 }
316 new_triplet = log_dens;
317
318 if(isnan(log_dens) || isinf(FABS(log_dens)))
319 {
320 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
321 MCMC_Print_Param(tree->mcmc,tree);
322 Exit("\n");
323 }
324 }
325 else
326 {
327 mu0 = mu1 = mu2 = -1.;
328 n0 = n1 = n2 = -1;
329
330 mu0 = tree->rates->br_r[n->num];
331 dt0 = FABS(tree->times->nd_t[n->num] - tree->times->nd_t[n->anc->num]);
332 n0 = tree->times->n_jps[n->num];
333 r0 = tree->rates->nd_r[n->num];
334
335 v1 = v2 = NULL;
336 for(i=0;i<3;i++)
337 {
338 if((n->v[i] != n->anc) && (n->b[i] != tree->e_root))
339 {
340 if(!v1)
341 {
342 v1 = n->v[i];
343 mu1 = tree->rates->br_r[v1->num];
344 dt1 = FABS(tree->times->nd_t[v1->num] - tree->times->nd_t[n->num]);
345 n1 = tree->times->n_jps[v1->num];
346 r1 = tree->rates->nd_r[v1->num];
347 }
348 else
349 {
350 v2 = n->v[i];
351 mu2 = tree->rates->br_r[v2->num];
352 dt2 = FABS(tree->times->nd_t[v2->num] - tree->times->nd_t[n->num]);
353 n2 = tree->times->n_jps[v2->num];
354 r2 = tree->rates->nd_r[v2->num];
355 }
356 }
357 }
358
359 mu1_mu0 = RATES_Lk_Rates_Core(mu0,mu1,r0,r1,n0,n1,dt0,dt1,tree);
360 mu2_mu0 = RATES_Lk_Rates_Core(mu0,mu2,r0,r1,n0,n2,dt0,dt2,tree);
361
362 new_triplet = mu1_mu0 + mu2_mu0;
363 }
364
365 tree->rates->c_lnL_rates = tree->rates->c_lnL_rates + new_triplet - curr_triplet;
366 tree->rates->triplet[n->num] = new_triplet;
367 }
368
369 //////////////////////////////////////////////////////////////
370 //////////////////////////////////////////////////////////////
371
372
373
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)374 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)
375 {
376 if((n1 > -1) && (n2 > -1))
377 {
378 return RATES_Compound_Core_Joint(mu1,mu2,n1,n2,dt1,dt2,alpha,beta,lexp,eps,approx);
379 }
380 else
381 {
382 return RATES_Compound_Core_Marginal(mu1,mu2,dt1,dt2,alpha,beta,lexp,eps,approx);
383 }
384 }
385
386 //////////////////////////////////////////////////////////////
387 //////////////////////////////////////////////////////////////
388
389
RATES_Compound_Core_Marginal(phydbl mu1,phydbl mu2,phydbl dt1,phydbl dt2,phydbl alpha,phydbl beta,phydbl lexp,phydbl eps,int approx)390 phydbl RATES_Compound_Core_Marginal(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, phydbl alpha, phydbl beta, phydbl lexp, phydbl eps, int approx)
391 {
392 phydbl p0,p1,v0,v1,v2;
393 phydbl dmu1;
394 int equ;
395 phydbl dens;
396
397 v0 = v1 = v2 = 0.0;
398
399 /* Probability of 0 and 1 jumps */
400 p0 = Dpois(0,lexp*(dt2+dt1),NO);
401 p1 = Dpois(1,lexp*(dt2+dt1),NO);
402
403 dmu1 = RATES_Dmu(mu1,-1,dt1,alpha,beta,lexp,0,0);
404
405 /* Are the two rates equal ? */
406 equ = 0;
407 if(FABS(mu1-mu2) < eps) equ = 1;
408
409 /* No jump */
410 if(equ)
411 {
412 v0 = 1.0*Dgamma(mu1,alpha,beta)/dmu1;
413 /* Rprintf("\n. mu1=%f mu2=%f",mu1,mu2); */
414 }
415 else
416 {
417 v0 = 1.E-100;
418 }
419
420 /* One jump */
421 v1 = RATES_Dmu1_And_Mu2_One_Jump_Two_Intervals(mu1,mu2,dt1,dt2,alpha,beta);
422 v1 /= dmu1;
423
424 /* Two jumps and more (approximation) */
425 if(approx == 1)
426 {
427 v2 =
428 RATES_Dmu(mu1,-1,dt1,alpha,beta,lexp,0,0)*RATES_Dmu(mu2,-1,dt2,alpha,beta,lexp,0,0) -
429 Dpois(0,lexp*dt1,NO) * Dpois(0,lexp*dt2,NO) *
430 Dgamma(mu1,alpha,beta) * Dgamma(mu2,alpha,beta);
431 }
432 else
433 {
434 v2 =
435 RATES_Dmu_One(mu1,dt1,alpha,beta,lexp) *
436 RATES_Dmu_One(mu2,dt2,alpha,beta,lexp);
437
438 v2 += Dpois(0,lexp*dt1,NO)*Dgamma(mu1,alpha,beta)*RATES_Dmu(mu2,-1,dt2,alpha,beta,lexp,1,0);
439 v2 += Dpois(0,lexp*dt2,NO)*Dgamma(mu2,alpha,beta)*RATES_Dmu(mu1,-1,dt1,alpha,beta,lexp,1,0);
440
441 }
442 /* PhyML_Printf("\n. %f %f %f %f %f ",mu1,mu2,dt1,dt2,v2); */
443 v2 /= dmu1;
444
445 dens = p0*v0 + p1*v1 + v2;
446 /* dens = p1*v1 + v2; */
447 /* dens = p1*v1 + v2; */
448 /* dens = v0; */
449 /* dens *= dmu1; */
450
451 if(dens < SMALL)
452 {
453 PhyML_Printf("\n. dens=%12G mu1=%12G mu2=%12G dt1=%12G dt2=%12G lexp=%12G alpha=%f v0=%f v1=%f v2=%f p0=%f p1=%f p2=%f",
454 dens,
455 mu1,mu2,dt1,dt2,
456 lexp,
457 alpha,
458 v0,v1,v2,
459 p0,p1,1.-p0-p1);
460 }
461
462 return dens;
463
464 }
465
466 /**********************************************************/
467
RATES_Compound_Core_Joint(phydbl mu1,phydbl mu2,int n1,int n2,phydbl dt1,phydbl dt2,phydbl alpha,phydbl beta,phydbl lexp,phydbl eps,int approx)468 phydbl RATES_Compound_Core_Joint(phydbl mu1, phydbl mu2, int n1, int n2, phydbl dt1, phydbl dt2,
469 phydbl alpha, phydbl beta, phydbl lexp, phydbl eps, int approx)
470 {
471 phydbl density;
472 phydbl dmu1;
473
474
475 if(n1 < 0 || n2 < 0)
476 {
477 PhyML_Fprintf(stderr,"\n. n1=%d n2=%d",n1,n2);
478 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
479 Warn_And_Exit("");
480 }
481
482 dmu1 = RATES_Dmu(mu1,n1,dt1,alpha,beta,lexp,0,0);
483
484 if((n1 < 0) || (n2 < 0))
485 {
486 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
487 Warn_And_Exit("");
488 }
489
490 if((n1 == 0) && (n2 == 0))
491 {
492 if(FABS(mu1-mu2) < eps) { density = Dgamma(mu1,alpha,beta); }
493 else { density = 1.E-70; }
494 }
495 else if((n1 == 0) && (n2 == 1))
496 {
497 density =
498 Dgamma(mu1,alpha,beta) *
499 RATES_Dmu1_Given_V_And_N(mu2,mu1,1,dt2,alpha,beta);
500 }
501 else if((n1 == 1) && (n2 == 0))
502 {
503 density =
504 Dgamma(mu2,alpha,beta) *
505 RATES_Dmu1_Given_V_And_N(mu1,mu2,1,dt1,alpha,beta);
506 }
507 else /* independent */
508 {
509 density =
510 RATES_Dmu(mu1,n1,dt1,alpha,beta,lexp,0,0) *
511 RATES_Dmu(mu2,n2,dt2,alpha,beta,lexp,0,0);
512 }
513
514 density /= dmu1;
515
516 density *= Dpois(n2,dt2*lexp,NO);
517
518 if(density < 1.E-70) density = 1.E-70;
519
520 /* PhyML_Printf("\n. density = %15G mu1=%3.4f mu2=%3.4f dt1=%3.4f dt2=%3.4f n1=%2d n2=%2d",density,mu1,mu2,dt1,dt2,n1,n2); */
521 return density;
522 }
523
524 /**********************************************************/
525
RATES_Print_Triplets(t_tree * tree)526 void RATES_Print_Triplets(t_tree *tree)
527 {
528 int i;
529 For(i,2*tree->n_otu-1) PhyML_Printf("\n. Node %3d t=%f",i,tree->rates->triplet[i]);
530 }
531
532
533 /**********************************************************/
534
RATES_Print_Rates(t_tree * tree)535 void RATES_Print_Rates(t_tree *tree)
536 {
537 RATES_Print_Rates_Pre(tree->n_root,tree->n_root->v[2],NULL,tree);
538 RATES_Print_Rates_Pre(tree->n_root,tree->n_root->v[1],NULL,tree);
539 }
540
541 //////////////////////////////////////////////////////////////
542 //////////////////////////////////////////////////////////////
543
RATES_Copy_Rate_Struct(t_rate * from,t_rate * to,int n_otu)544 void RATES_Copy_Rate_Struct(t_rate *from, t_rate *to, int n_otu)
545 {
546 int i;
547
548 to->lexp = from->lexp;
549 to->alpha = from->alpha;
550 to->less_likely = from->less_likely;
551
552 to->nu = from->nu;
553 to->min_nu = from->min_nu;
554 to->max_nu = from->max_nu;
555
556 to->min_rate = from->min_rate;
557 to->max_rate = from->max_rate;
558
559 to->c_lnL1 = from->c_lnL1;
560 to->c_lnL2 = from->c_lnL2;
561
562 to->c_lnL_rates = from->c_lnL_rates;
563
564 to->clock_r = from->clock_r;
565 to->min_clock = from->min_clock;
566 to->max_clock = from->max_clock;
567
568 to->lbda_nu = from->lbda_nu;
569 to->min_dt = from->min_dt;
570 to->step_rate = from->step_rate;
571 to->true_tree_size = from->true_tree_size;
572 to->p_max = from->p_max;
573 to->norm_fact = from->norm_fact;
574
575 to->adjust_rates = from->adjust_rates;
576 to->use_rates = from->use_rates;
577 to->bl_from_rt = from->bl_from_rt;
578 to->approx = from->approx;
579 to->model = from->model;
580 to->is_allocated = from->is_allocated;
581 to->met_within_gibbs = from->met_within_gibbs;
582
583 to->update_mean_l = from->update_mean_l;
584 to->update_cov_l = from->update_cov_l;
585
586 to->br_r_recorded = from->br_r_recorded;
587
588 to->log_K_cur = from->log_K_cur;
589 to->cur_comb_numb = from->cur_comb_numb;
590
591 for(i=0;i<2*n_otu-1;++i) to->nd_r[i] = from->nd_r[i];
592 for(i=0;i<2*n_otu-1;++i) to->br_r[i] = from->br_r[i];
593 for(i=0;i<2*n_otu-1;++i) to->buff_br_r[i] = from->buff_br_r[i];
594 for(i=0;i<2*n_otu-1;++i) to->buff_nd_r[i] = from->buff_nd_r[i];
595 for(i=0;i<2*n_otu-1;++i) to->true_r[i] = from->true_r[i];
596 for(i=0;i<2*n_otu-2;++i) to->dens[i] = from->dens[i];
597 for(i=0;i<2*n_otu-1;++i) to->triplet[i] = from->triplet[i];
598 For(i,(2*n_otu-2)*(2*n_otu-2)) to->cov_l[i] = from->cov_l[i];
599 For(i,(2*n_otu-2)*(2*n_otu-2)) to->invcov[i] = from->invcov[i];
600 for(i=0;i<2*n_otu-2;++i) to->mean_l[i] = from->mean_l[i];
601 for(i=0;i<2*n_otu-2;++i) to->ml_l[i] = from->ml_l[i];
602 for(i=0;i<2*n_otu-2;++i) to->cur_l[i] = from->cur_l[i];
603 For(i,2*n_otu-3) to->u_ml_l[i] = from->u_ml_l[i];
604 For(i,2*n_otu-3) to->u_cur_l[i] = from->u_cur_l[i];
605 For(i,(2*n_otu-2)*(2*n_otu-2)) to->cov_r[i] = from->cov_r[i];
606 for(i=0;i<2*n_otu-2;++i) to->cond_var[i] = from->cond_var[i];
607 for(i=0;i<2*n_otu-2;++i) to->mean_r[i] = from->mean_r[i];
608 For(i,(2*n_otu-1)*(2*n_otu-1)) to->lca[i] = from->lca[i];
609 For(i,(2*n_otu-2)*(2*n_otu-2)) to->reg_coeff[i] = from->reg_coeff[i];
610 For(i,(2*n_otu-2)*(6*n_otu-9)) to->trip_reg_coeff[i] = from->trip_reg_coeff[i];
611 For(i,(2*n_otu-2)*9) to->trip_cond_cov[i] = from->trip_cond_cov[i];
612 For(i,2*n_otu) to->_2n_vect1[i] = from->_2n_vect1[i];
613 For(i,2*n_otu) to->_2n_vect2[i] = from->_2n_vect2[i];
614 For(i,2*n_otu) to->_2n_vect3[i] = from->_2n_vect3[i];
615 For(i,2*n_otu) to->_2n_vect4[i] = from->_2n_vect4[i];
616 For(i,2*n_otu) to->_2n_vect5[i] = from->_2n_vect5[i];
617 For(i,4*n_otu*n_otu) to->_2n2n_vect1[i] = from->_2n2n_vect1[i];
618 For(i,4*n_otu*n_otu) to->_2n2n_vect2[i] = from->_2n2n_vect2[i];
619 for(i=0;i<2*n_otu-1;++i) to->br_do_updt[i] = from->br_do_updt[i];
620 for(i=0;i<2*n_otu-1;++i) to->cur_gamma_prior_mean[i] = from->cur_gamma_prior_mean[i];
621 for(i=0;i<2*n_otu-1;++i) to->cur_gamma_prior_var[i] = from->cur_gamma_prior_var[i];
622 for(i=0;i<2*n_otu-1;++i) to->n_tips_below[i] = from->n_tips_below[i];
623 }
624
625 //////////////////////////////////////////////////////////////
626 //////////////////////////////////////////////////////////////
627
628
RATES_Duplicate_Calib_Struct(t_tree * from,t_tree * to)629 void RATES_Duplicate_Calib_Struct(t_tree *from, t_tree *to)
630 {
631 int i,j;
632
633 to->times->n_cal = from->times->n_cal;
634 for(i=0;i<from->times->n_cal;i++)
635 {
636 to->times->a_cal[i] = Duplicate_Calib(from->times->a_cal[i]);
637 for(j=0;j<from->times->a_cal[i]->clade_list_size;j++)
638 Init_Target_Tip(to->times->a_cal[i]->clade_list[j],to);
639 }
640 }
641
642 //////////////////////////////////////////////////////////////
643 //////////////////////////////////////////////////////////////
644
RATES_Print_Rates_Pre(t_node * a,t_node * d,t_edge * b,t_tree * tree)645 void RATES_Print_Rates_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
646 {
647 if((d == tree->n_root->v[2] && d->tax) || (d == tree->n_root->v[1] && d->tax))
648 PhyML_Printf("\n. a=%3d ++d=%3d rate=%12f t_left=%12f t_rght=%12f ml=%12f l=%12f %12f",
649 a->num,d->num,
650 tree->rates->br_r[d->num],
651 tree->times->nd_t[a->num],tree->times->nd_t[d->num],
652 tree->rates->ml_l[d->num],
653 tree->rates->cur_l[d->num],
654 (tree->times->nd_t[d->num]-tree->times->nd_t[a->num])*tree->rates->clock_r*tree->rates->br_r[d->num]);
655
656 else if((d == tree->n_root->v[2]) || (d == tree->n_root->v[1]))
657 PhyML_Printf("\n. a=%3d __d=%3d rate=%12f t_left=%12f t_rght=%12f ml=%12f l=%12f %12f",
658 a->num,d->num,
659 tree->rates->br_r[d->num],
660 tree->times->nd_t[a->num],tree->times->nd_t[d->num],
661 tree->rates->ml_l[d->num],
662 tree->rates->cur_l[d->num],
663 (tree->times->nd_t[d->num]-tree->times->nd_t[a->num])*tree->rates->clock_r*tree->rates->br_r[d->num]);
664 else
665 PhyML_Printf("\n. a=%3d d=%3d rate=%12f t_left=%12f t_rght=%12f ml=%12f l=%12f %12f",
666 a->num,d->num,
667 tree->rates->br_r[d->num],
668 tree->times->nd_t[a->num],tree->times->nd_t[d->num],
669 tree->rates->ml_l[d->num],
670 tree->rates->cur_l[d->num],
671 (tree->times->nd_t[d->num]-tree->times->nd_t[a->num])*tree->rates->clock_r*tree->rates->br_r[d->num]);
672
673 if(d->tax) return;
674 else
675 {
676 int i;
677
678 for(i=0;i<3;i++)
679 {
680 if((d->v[i] != a) && (d->b[i] != tree->e_root))
681 {
682 RATES_Print_Rates_Pre(d,d->v[i],d->b[i],tree);
683 }
684 }
685 }
686 }
687
688 //////////////////////////////////////////////////////////////
689 //////////////////////////////////////////////////////////////
690
691
RATES_Average_Rate(t_tree * tree)692 phydbl RATES_Average_Rate(t_tree *tree)
693 {
694 int i;
695 phydbl sum;
696 sum = 0.0;
697 for(i=0;i<2*tree->n_otu-2;++i) sum += tree->rates->br_r[i];
698 return sum/(2*tree->n_otu-2);
699 }
700 //////////////////////////////////////////////////////////////
701 //////////////////////////////////////////////////////////////
702
RATES_Average_Substitution_Rate(t_tree * tree)703 phydbl RATES_Average_Substitution_Rate(t_tree *tree)
704 {
705 phydbl l,dt;
706 int i;
707
708 l = 0.0;
709 dt = 0.0;
710 for(i=0;i<2*tree->n_otu-1;++i)
711 {
712 if(tree->a_nodes[i] != tree->n_root)
713 {
714 dt += FABS(tree->times->nd_t[tree->a_nodes[i]->num] - tree->times->nd_t[tree->a_nodes[i]->anc->num]);
715 l += tree->rates->cur_l[tree->a_nodes[i]->num];
716 }
717 }
718
719 return(l/dt);
720 }
721
722 //////////////////////////////////////////////////////////////
723 //////////////////////////////////////////////////////////////
724
RATES_Check_Mean_Rates_True(t_tree * tree)725 phydbl RATES_Check_Mean_Rates_True(t_tree *tree)
726 {
727 phydbl sum;
728 int i;
729
730 sum = 0.0;
731 for(i=0;i<2*tree->n_otu-2;++i) sum += tree->rates->true_r[i];
732 return(sum/(phydbl)(2*tree->n_otu-2));
733 }
734
735 //////////////////////////////////////////////////////////////
736 //////////////////////////////////////////////////////////////
737
RATES_Check_Node_Times(t_tree * tree)738 int RATES_Check_Node_Times(t_tree *tree)
739 {
740 int err;
741 err = NO;
742 RATES_Check_Node_Times_Pre(tree->n_root,tree->n_root->v[2],&err,tree);
743 RATES_Check_Node_Times_Pre(tree->n_root,tree->n_root->v[1],&err,tree);
744 return err;
745 }
746
747 //////////////////////////////////////////////////////////////
748 //////////////////////////////////////////////////////////////
749
RATES_Check_Node_Times_Pre(t_node * a,t_node * d,int * err,t_tree * tree)750 void RATES_Check_Node_Times_Pre(t_node *a, t_node *d, int *err, t_tree *tree)
751 {
752 if((tree->times->nd_t[d->num] < tree->times->nd_t[a->num]) || (FABS(tree->times->nd_t[d->num] - tree->times->nd_t[a->num]) < 1.E-20))
753 {
754 PhyML_Printf("\n. a->t=%f d->t=%f",tree->times->nd_t[a->num],tree->times->nd_t[d->num]);
755 PhyML_Printf("\n. a->t_prior_min=%f a->t_prior_max=%f",tree->times->t_prior_min[a->num],tree->times->t_prior_max[a->num]);
756 PhyML_Printf("\n. d->t_prior_min=%f d->t_prior_max=%f",tree->times->t_prior_min[d->num],tree->times->t_prior_max[d->num]);
757 *err = YES;
758 }
759 if(d->tax) return;
760 else
761 {
762 int i;
763
764 for(i=0;i<3;i++)
765 if((d->v[i] != a) && (d->b[i] != tree->e_root))
766 RATES_Check_Node_Times_Pre(d,d->v[i],err,tree);
767 }
768 }
769 //////////////////////////////////////////////////////////////
770 //////////////////////////////////////////////////////////////
771
RATES_Bracket_N_Jumps(int * up,int * down,phydbl param)772 void RATES_Bracket_N_Jumps(int *up, int *down, phydbl param)
773 {
774 phydbl cdf,eps,a,b,c;
775 int step;
776
777 step = 10;
778 eps = 1.E-10;
779 cdf = 0.0;
780 c = 1;
781
782 while(cdf < 1.-eps)
783 {
784 c = (int)FLOOR(c * step);
785 cdf = Ppois(c,param);
786 }
787
788 a = 0.0;
789 b = (c-a)/2.;
790 step = 0;
791 do
792 {
793 step++;
794 cdf = Ppois(b,param);
795 if(cdf < eps) a = b;
796 else
797 {
798 break;
799 }
800 b = (c-a)/2.;
801 }
802 while(step < 1000);
803
804 if(step == 1000)
805 {
806 PhyML_Fprintf(stderr,"\n. a=%f b=%f c=%f param=%f",a,b,c,param);
807 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
808 Warn_And_Exit("");
809 }
810 *up = c;
811 *down = a;
812 }
813
814 //////////////////////////////////////////////////////////////
815 //////////////////////////////////////////////////////////////
816
817 /*
818 mu : average rate of the time period dt
819 dt : time period to be considered
820 a : rate at a given time point is gamma distributed. a is the shape parameter
821 b : rate at a given time point is gamma distributed. b is the scale parameter
822 lexp : the number of rate switches is Poisson distributed with parameter lexp * dt
823 */
824 /* compute f(mu;dt,a,b,lexp), the probability density of mu. We need to integrate over the
825 possible number of jumps (n) during the time interval dt */
RATES_Dmu(phydbl mu,int n_jumps,phydbl dt,phydbl a,phydbl b,phydbl lexp,int min_n,int jps_dens)826 phydbl RATES_Dmu(phydbl mu, int n_jumps, phydbl dt, phydbl a, phydbl b, phydbl lexp, int min_n, int jps_dens)
827 {
828 if(n_jumps < 0) /* Marginal, i.e., the number of jumps is not fixed */
829 {
830 phydbl var,cumpoissprob,dens,mean,poissprob,ab2,gammadens,lexpdt;
831 int n,up,down;
832
833 var = 0.0;
834 cumpoissprob = 0.0;
835 dens = 0.0;
836 n = 0;
837 mean = a*b;
838 ab2 = a*b*b;
839 lexpdt = lexp*dt;
840
841 RATES_Bracket_N_Jumps(&up,&down,lexpdt);
842 For(n,MAX(down,min_n)-1) cumpoissprob += Dpois(n,lexpdt,NO);
843
844 for(n=MAX(down,min_n);n<up+1;n++)
845 {
846 poissprob = Dpois(n,lexpdt,NO); /* probability of having n jumps */
847 var = (2./(n+2.))*ab2; /* var(mu|n) = var(mu|n=0) * 2 / (n+2) */
848 gammadens = Dgamma_Moments(mu,mean,var);
849 dens += poissprob * gammadens;
850 cumpoissprob += poissprob;
851 if(cumpoissprob > 1.-1.E-04) break;
852 }
853
854 if(dens < 1.E-70) dens = 1.E-70;
855
856 return(dens);
857 }
858 else /* Joint, i.e., return P(mu | dt, n_jumps) */
859 {
860 phydbl mean, var, density;
861
862
863 mean = 1.0;
864 var = (2./(n_jumps+2.))*a*b*b;
865
866 if(jps_dens)
867 density = Dgamma_Moments(mu,mean,var) * Dpois(n_jumps,dt*lexp,NO);
868 else
869 density = Dgamma_Moments(mu,mean,var);
870
871 if(density < 1.E-70) density = 1.E-70;
872
873 return density;
874 }
875 }
876
877 //////////////////////////////////////////////////////////////
878 //////////////////////////////////////////////////////////////
879
880
RATES_Dmu_One(phydbl mu,phydbl dt,phydbl a,phydbl b,phydbl lexp)881 phydbl RATES_Dmu_One(phydbl mu, phydbl dt, phydbl a, phydbl b, phydbl lexp)
882 {
883 phydbl var,cumpoissprob,dens,mean,poissprob,ab2,gammadens,lexpdt;
884 int n,up,down;
885
886 var = 0.0;
887 cumpoissprob = 0.0;
888 dens = 0.0;
889 n = 0;
890 mean = a*b;
891 ab2 = a*b*b;
892 lexpdt = lexp*dt;
893
894 if(dt < 0.0)
895 {
896 PhyML_Fprintf(stderr,"\n. dt=%f",dt);
897 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
898 Warn_And_Exit("");
899 }
900
901 if(lexpdt < SMALL)
902 {
903 PhyML_Fprintf(stderr,"\n. lexpdt=%G lexp=%G dt=%G",lexpdt,lexp,dt);
904 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
905 Warn_And_Exit("");
906 }
907
908 if(mu < 1.E-10)
909 {
910 PhyML_Fprintf(stderr,"\n. mu=%G",mu);
911 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
912 Warn_And_Exit("");
913 }
914
915 RATES_Bracket_N_Jumps(&up,&down,lexpdt);
916
917 For(n,MAX(1,down)-1) cumpoissprob += Dpois(n,lexpdt,NO);
918
919 for(n=MAX(1,down);n<up+1;n++) /* WARNING: we are considering that at least one jump occurs in the interval */
920 {
921 poissprob = Dpois(n,lexpdt,NO); /* probability of having n jumps */
922 var = (n/((n+1)*(n+1)*(n+2)))*(POW(1-a*b,2) + 2/(n+1)*ab2) + 2*n*n*ab2/POW(n+1,3);
923 gammadens = Dgamma_Moments(mu,mean,var);
924 dens += poissprob * gammadens;
925 cumpoissprob += poissprob;
926 if(cumpoissprob > 1.-1.E-06) break;
927 }
928
929 return(dens);
930 }
931
932 //////////////////////////////////////////////////////////////
933 //////////////////////////////////////////////////////////////
934
935 /* Given the times of nodes a (ta) and d (td), the shape of the gamma distribution of instantaneous
936 rates, the parameter of the exponential distribution of waiting times between rate jumps and the
937 instantaneous rate at t_node a, this function works out an expected number of (amino-acids or
938 nucleotide) substitutions per site.
939 */
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)940 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)
941 {
942 phydbl curr_r, curr_t, next_t;
943
944 switch(rates->model)
945 {
946 case COMPOUND_COR:case COMPOUND_NOCOR:
947 {
948 /* Compound Poisson */
949 if(rates->model == COMPOUND_COR)
950 {
951 curr_r = r_beg;
952 *mean_r = r_beg;
953 }
954 else
955 {
956 curr_r = Rgamma(rates->nu,1./rates->nu);;
957 *mean_r = curr_r;
958 }
959
960 curr_t = t_beg + Rexp(rates->lexp); /* Exponentially distributed waiting times */
961 next_t = curr_t;
962
963 *n_jumps = 0;
964 while(curr_t < t_end)
965 {
966 curr_r = Rgamma(rates->nu,1./rates->nu); /* Gamma distributed random instantaneous rate */
967
968 (*n_jumps)++;
969
970 next_t = curr_t + Rexp(rates->lexp);
971
972 if(next_t < t_end)
973 {
974 *mean_r = (1./(next_t - t_beg)) * (*mean_r * (curr_t - t_beg) + curr_r * (next_t - curr_t));
975 }
976 else
977 {
978 *mean_r = (1./(t_end - t_beg)) * (*mean_r * (curr_t - t_beg) + curr_r * (t_end - curr_t));
979 }
980 curr_t = next_t;
981 }
982
983 /* PhyML_Printf("\n. [%3d %f %f]",*n_jumps,*mean_r,r_beg); */
984
985 if(*mean_r < rates->min_rate) *mean_r = rates->min_rate;
986 if(*mean_r > rates->max_rate) *mean_r = rates->max_rate;
987
988 *r_end = curr_r;
989 break;
990 }
991 case EXPONENTIAL:
992 {
993 *mean_r = Rexp(rates->nu);
994
995 if(*mean_r < rates->min_rate) *mean_r = rates->min_rate;
996 if(*mean_r > rates->max_rate) *mean_r = rates->max_rate;
997
998 *r_end = *mean_r;
999 break;
1000 }
1001 case LOGNORMAL:
1002 {
1003 *mean_r = Rgamma(rates->nu,1./rates->nu);
1004
1005 if(*mean_r < rates->min_rate) *mean_r = rates->min_rate;
1006 if(*mean_r > rates->max_rate) *mean_r = rates->max_rate;
1007
1008 *r_end = *mean_r;
1009 break;
1010 }
1011 case THORNE:
1012 {
1013 phydbl sd,mean;
1014 int err;
1015
1016 sd = SQRT(rates->nu*FABS(t_beg-t_end));
1017 mean = r_beg;
1018
1019 *mean_r = Rnorm_Trunc(mean,sd,rates->min_rate,rates->max_rate,&err);
1020
1021 if(err) PhyML_Printf("\n. %s %d %d",__FILE__,__LINE__,tree->mcmc->run);
1022 *r_end = *mean_r;
1023 break;
1024 }
1025 default:
1026 {
1027 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
1028 Exit("\n. Model not implemented yet.\n");
1029 break;
1030 }
1031 }
1032 }
1033
1034 //////////////////////////////////////////////////////////////
1035 //////////////////////////////////////////////////////////////
1036
1037
RATES_Get_Mean_Rates_Pre(t_node * a,t_node * d,t_edge * b,phydbl r_a,t_tree * tree)1038 void RATES_Get_Mean_Rates_Pre(t_node *a, t_node *d, t_edge *b, phydbl r_a, t_tree *tree)
1039 {
1040 phydbl a_t,d_t;
1041 phydbl mean_r;
1042 int n_jumps;
1043 phydbl r_d;
1044
1045 a_t = tree->times->nd_t[a->num];
1046 d_t = tree->times->nd_t[d->num];
1047
1048 mean_r = -1.;
1049 n_jumps = -1;
1050 r_d = -1.;
1051 RATES_Expect_Number_Subst(a_t,d_t,r_a,&n_jumps,&mean_r,&r_d,tree->rates,tree);
1052
1053 tree->rates->br_r[d->num] = mean_r;
1054 tree->rates->true_r[d->num] = mean_r;
1055 tree->times->t_jps[d->num] = n_jumps;
1056
1057
1058 /* Move to the next branches */
1059 if(d->tax) return;
1060 else
1061 {
1062 int i;
1063
1064 for(i=0;i<3;i++)
1065 {
1066 if((d->v[i] != a) && (d->b[i] != tree->e_root))
1067 {
1068 RATES_Get_Mean_Rates_Pre(d,d->v[i],d->b[i],r_d,tree);
1069 }
1070 }
1071 }
1072
1073 }
1074
1075 //////////////////////////////////////////////////////////////
1076 //////////////////////////////////////////////////////////////
1077
RATES_Random_Branch_Lengths(t_tree * tree)1078 void RATES_Random_Branch_Lengths(t_tree *tree)
1079 {
1080 phydbl r0;
1081
1082 r0 = 1.0;
1083
1084 tree->rates->br_r[tree->n_root->num] = r0;
1085
1086 RATES_Get_Mean_Rates_Pre(tree->n_root,tree->n_root->v[2],NULL,r0,tree);
1087 RATES_Get_Mean_Rates_Pre(tree->n_root,tree->n_root->v[1],NULL,r0,tree);
1088
1089 /* RATES_Normalise_Rates(tree); */
1090
1091 RATES_Update_Cur_Bl(tree);
1092 RATES_Initialize_True_Rates(tree);
1093
1094 tree->n_root_pos =
1095 tree->rates->cur_l[tree->n_root->v[2]->num] /
1096 (tree->rates->cur_l[tree->n_root->v[2]->num] + tree->rates->cur_l[tree->n_root->v[1]->num]);
1097 }
1098
1099 //////////////////////////////////////////////////////////////
1100 //////////////////////////////////////////////////////////////
1101
1102
RATES_Set_Node_Times(t_tree * tree)1103 void RATES_Set_Node_Times(t_tree *tree)
1104 {
1105 RATES_Set_Node_Times_Pre(tree->n_root,tree->n_root->v[2],tree);
1106 RATES_Set_Node_Times_Pre(tree->n_root,tree->n_root->v[1],tree);
1107 }
1108
1109 //////////////////////////////////////////////////////////////
1110 //////////////////////////////////////////////////////////////
1111
1112
RATES_Init_Triplets(t_tree * tree)1113 void RATES_Init_Triplets(t_tree *tree)
1114 {
1115 int i;
1116 For(i,2*tree->n_otu-1) tree->rates->triplet[i] = 0.0;
1117 }
1118 //////////////////////////////////////////////////////////////
1119 //////////////////////////////////////////////////////////////
1120
1121
RATES_Set_Node_Times_Pre(t_node * a,t_node * d,t_tree * tree)1122 void RATES_Set_Node_Times_Pre(t_node *a, t_node *d, t_tree *tree)
1123 {
1124 if(d->tax) return;
1125 else
1126 {
1127 t_node *v1, *v2; /* the two sons of d */
1128 phydbl t_sup, t_inf;
1129 int i;
1130
1131 v1 = v2 = NULL;
1132 for(i=0;i<3;i++) if((d->v[i] != a) && (d->b[i] != tree->e_root))
1133 {
1134 if(!v1) v1 = d->v[i];
1135 else v2 = d->v[i];
1136 }
1137
1138 t_inf = MIN(tree->times->nd_t[v1->num],tree->times->nd_t[v2->num]);
1139 t_sup = tree->times->nd_t[a->num];
1140
1141 if(t_sup > t_inf)
1142 {
1143 PhyML_Fprintf(stderr,"\n. t_sup = %f t_inf = %f",t_sup,t_inf);
1144 PhyML_Fprintf(stderr,"\n. Run = %d",tree->mcmc->run);
1145 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1146 Warn_And_Exit("");
1147 }
1148
1149 if(tree->times->nd_t[d->num] > t_inf)
1150 {
1151 tree->times->nd_t[d->num] = t_inf;
1152 PhyML_Printf("\n. Time for t_node %d is larger than should be. Set it to %f",d->num,tree->times->nd_t[d->num]);
1153 }
1154 else if(tree->times->nd_t[d->num] < t_sup)
1155 {
1156 tree->times->nd_t[d->num] = t_sup;
1157 PhyML_Printf("\n. Time for t_node %d is lower than should be. Set it to %f",d->num,tree->times->nd_t[d->num]);
1158 }
1159
1160 for(i=0;i<3;i++)
1161 {
1162 if((d->v[i] != a) && (d->b[i] != tree->e_root))
1163 {
1164 RATES_Set_Node_Times_Pre(d,d->v[i],tree);
1165 }
1166 }
1167 }
1168 }
1169
1170 //////////////////////////////////////////////////////////////
1171 //////////////////////////////////////////////////////////////
1172
1173
1174
RATES_Dmu1_And_Mu2_One_Jump_Two_Intervals(phydbl mu1,phydbl mu2,phydbl dt1,phydbl dt2,phydbl alpha,phydbl beta)1175 phydbl RATES_Dmu1_And_Mu2_One_Jump_Two_Intervals(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, phydbl alpha, phydbl beta)
1176 {
1177 phydbl dens;
1178
1179 if(mu2 < 1.E-10)
1180 {
1181 PhyML_Fprintf(stderr,"\n. mu2=%G",mu2);
1182 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
1183 Warn_And_Exit("");
1184 }
1185
1186 if(mu1 < 1.E-10)
1187 {
1188 PhyML_Fprintf(stderr,"\n. mu2=%G",mu1);
1189 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1190 Warn_And_Exit("");
1191 }
1192
1193 dens =
1194 ((dt1/(dt1+dt2)) * RATES_Dmu1_Given_V_And_N(mu1,mu2,1,dt1,alpha,beta) * Dgamma(mu2,alpha,beta)) +
1195 ((dt2/(dt1+dt2)) * RATES_Dmu1_Given_V_And_N(mu2,mu1,1,dt2,alpha,beta) * Dgamma(mu1,alpha,beta));
1196
1197 return dens;
1198 }
1199
1200 //////////////////////////////////////////////////////////////
1201 //////////////////////////////////////////////////////////////
1202
1203
RATES_Dmu1_Given_V_And_N(phydbl mu1,phydbl v,int n,phydbl dt1,phydbl a,phydbl b)1204 phydbl RATES_Dmu1_Given_V_And_N(phydbl mu1, phydbl v, int n, phydbl dt1, phydbl a, phydbl b)
1205 {
1206 phydbl lbda,dens,h,u;
1207 phydbl mean,var;
1208 int n_points,i;
1209 phydbl ndb;
1210 phydbl end, beg;
1211
1212 n_points = 20;
1213
1214 end = MIN(mu1/v-0.01,0.99);
1215 beg = 0.01;
1216
1217 dens = 0.0;
1218
1219 if(end > beg)
1220 {
1221 mean = a*b;
1222 var = a*b*b*2./(n+1.);
1223 ndb = (phydbl)n/dt1;
1224
1225 h = (end - beg) / (phydbl)n_points;
1226
1227 lbda = beg;
1228 for(i=0;i<n_points-1;i++)
1229 {
1230 lbda += h;
1231 u = (mu1 - lbda*v)/(1.-lbda);
1232
1233 if(u < 1.E-10)
1234 {
1235 PhyML_Fprintf(stderr,"\n. u = %G",u);
1236 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1237 Warn_And_Exit("");
1238 }
1239
1240 dens += Dgamma_Moments(u,mean,var) / (1.-lbda) * ndb * POW(1.-lbda,n-1);
1241 }
1242 dens *= 2.;
1243
1244 lbda = beg;
1245 u = (mu1 - lbda*v)/(1.-lbda);
1246 if(u < 1.E-10)
1247 {
1248 PhyML_Fprintf(stderr,"\n. mu1 = %f lambda = %f v = %f u = %G beg = %f end = %f",mu1,lbda,v,u,beg,end);
1249 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1250 Warn_And_Exit("");
1251 }
1252
1253 dens += Dgamma_Moments(u,mean,var) / (1.-lbda) * ndb * POW(1.-lbda,n-1);
1254
1255 lbda = end;
1256 u = (mu1 - lbda*v)/(1.-lbda);
1257 if(u < 1.E-10)
1258 {
1259 PhyML_Fprintf(stderr,"\n. u = %G",u);
1260 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1261 Warn_And_Exit("");
1262 }
1263
1264 dens += Dgamma_Moments(u,mean,var) / (1.-lbda) * ndb * POW(1.-lbda,n-1);
1265
1266 dens *= (h/2.);
1267 dens *= dt1;
1268 }
1269
1270 return(dens);
1271 }
1272
1273 //////////////////////////////////////////////////////////////
1274 //////////////////////////////////////////////////////////////
1275
1276 /* Joint density of mu1 and a minimum number of jumps occuring in the interval dt1+dt2 given mu1.
1277 1 jump occurs at the junction of the two intervals, which makes mu1 and mu2 independant */
RATES_Dmu2_And_Mu1_Given_Min_N(phydbl mu1,phydbl mu2,phydbl dt1,phydbl dt2,int n_min,phydbl a,phydbl b,phydbl lexp)1278 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)
1279 {
1280 phydbl density, lexpdt,cumpoiss,poiss;
1281 int i;
1282 int up,down;
1283
1284 density = 0.0;
1285 lexpdt = lexp * (dt1+dt2);
1286 cumpoiss = 0.0;
1287
1288 RATES_Bracket_N_Jumps(&up,&down,lexpdt);
1289
1290 For(i,MAX(up,n_min)-1)
1291 {
1292 poiss = Dpois(i,lexpdt,NO);
1293 cumpoiss = cumpoiss + poiss;
1294 }
1295
1296 for(i=MAX(up,n_min);i<up;i++)
1297 {
1298 /* poiss = Dpois(i-1,lexpdt,NO); /\* Complies with the no correlation model *\/ */
1299 poiss = Dpois(i,lexpdt,NO);
1300 cumpoiss = cumpoiss + poiss;
1301
1302 density = density + poiss * RATES_Dmu2_And_Mu1_Given_N(mu1,mu2,dt1,dt2,i-1,a,b,lexp);
1303 /* density = density + poiss * RATES_Dmu2_And_Mu1_Given_N_Full(mu1,mu2,dt1,dt2,i,a,b,lexp); */
1304 if(cumpoiss > 1.-1.E-6) break;
1305 }
1306
1307 if(density < 0.0)
1308 {
1309 PhyML_Fprintf(stderr,"\n. density=%f cmpoiss = %f i=%d n_min=%d mu1=%f mu2=%f dt1=%f dt2=%f",
1310 density,cumpoiss,i,n_min,mu1,mu2,dt1,dt2);
1311 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1312 Warn_And_Exit("");
1313 }
1314
1315 return(density);
1316 }
1317
1318 //////////////////////////////////////////////////////////////
1319 //////////////////////////////////////////////////////////////
1320
1321 /* Joint density of mu1 and mu2 given the number of jumps (n) in the interval dt1+dt2, which are considered
1322 as independant. Hence, for n jumps occuring in dt1+dt2, the number of jumps occuring in dt1 and dt2 are
1323 n and 0, or n-1 and 1, or n-2 and 2 ... or 0 and n. This function sums over all these scenarios, with
1324 weights corresponding to the probability of each partitition.
1325 */
1326
RATES_Dmu2_And_Mu1_Given_N(phydbl mu1,phydbl mu2,phydbl dt1,phydbl dt2,int n,phydbl a,phydbl b,phydbl lexp)1327 phydbl RATES_Dmu2_And_Mu1_Given_N(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, int n, phydbl a, phydbl b, phydbl lexp)
1328 {
1329 phydbl density,cumpoiss,poiss,abb,ab,lognf,logdt1,logdt2,nlogdt;
1330 int i;
1331
1332 density = 0.0;
1333 abb = a*b*b;
1334 ab = a*b;
1335 cumpoiss = 0.0;
1336 poiss = 0.0;
1337 lognf = LnFact(n);
1338 logdt1 = log(dt1);
1339 logdt2 = log(dt2);
1340 nlogdt = n*log(dt1+dt2);
1341
1342 For(i,n+1)
1343 {
1344 poiss = lognf - LnFact(i) - LnFact(n-i) + i*logdt1 + (n-i)*logdt2 - nlogdt;
1345 poiss = exp(poiss);
1346 cumpoiss = cumpoiss + poiss;
1347
1348 if(mu2 < 1.E-10)
1349 {
1350 PhyML_Fprintf(stderr,"\n. mu2=%f",mu2);
1351 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1352 Warn_And_Exit("");
1353 }
1354
1355 if(mu1 < 1.E-10)
1356 {
1357 PhyML_Fprintf(stderr,"\n. mu1=%f",mu1);
1358 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1359 Warn_And_Exit("");
1360 }
1361
1362
1363 density = density + poiss * Dgamma_Moments(mu2,ab,2./((n-i)+2.)*abb) * Dgamma_Moments(mu1,ab,2./(i+2.)*abb);
1364 if(cumpoiss > 1.-1.E-6) break;
1365 }
1366
1367 return(density);
1368 }
1369
1370 //////////////////////////////////////////////////////////////
1371 //////////////////////////////////////////////////////////////
1372
RATES_Initialize_True_Rates(t_tree * tree)1373 void RATES_Initialize_True_Rates(t_tree *tree)
1374 {
1375 int i;
1376 For(i,2*tree->n_otu-2) tree->rates->true_r[i] = tree->rates->br_r[i];
1377 }
1378 //////////////////////////////////////////////////////////////
1379 //////////////////////////////////////////////////////////////
1380
1381
RATES_Get_Rates_From_Bl(t_tree * tree)1382 void RATES_Get_Rates_From_Bl(t_tree *tree)
1383 {
1384 phydbl dt,cr;
1385 t_node *left, *rght;
1386 int i;
1387
1388 dt = -1.0;
1389 cr = tree->rates->clock_r;
1390
1391 if(tree->n_root)
1392 {
1393 dt = FABS(tree->times->nd_t[tree->n_root->num] - tree->times->nd_t[tree->n_root->v[2]->num]);
1394 tree->rates->br_r[tree->n_root->v[2]->num] = 0.5 * tree->e_root->l->v / (dt*cr);
1395 dt = FABS(tree->times->nd_t[tree->n_root->num] - tree->times->nd_t[tree->n_root->v[1]->num]);
1396 tree->rates->br_r[tree->n_root->v[1]->num] = 0.5 * tree->e_root->l->v / (dt*cr);
1397 }
1398
1399
1400 For(i,2*tree->n_otu-3)
1401 {
1402 if(tree->a_edges[i] != tree->e_root)
1403 {
1404 left = tree->a_edges[i]->left;
1405 rght = tree->a_edges[i]->rght;
1406 dt = FABS(tree->times->nd_t[left->num] - tree->times->nd_t[rght->num]);
1407
1408 if(left->anc == rght) tree->rates->br_r[left->num] = tree->a_edges[i]->l->v / (dt*cr);
1409 else tree->rates->br_r[rght->num] = tree->a_edges[i]->l->v / (dt*cr);
1410 }
1411 }
1412 }
1413
1414 //////////////////////////////////////////////////////////////
1415 //////////////////////////////////////////////////////////////
1416
RATES_Lk_Jumps(t_tree * tree)1417 phydbl RATES_Lk_Jumps(t_tree *tree)
1418 {
1419 int i,n_jps;
1420 phydbl dens,dt,lexp;
1421 t_node *n;
1422
1423 n = NULL;
1424 lexp = tree->rates->lexp;
1425 n_jps = 0;
1426 dt = 0.0;
1427 dens = 0.0;
1428
1429 for(i=0;i<2*tree->n_otu-2;++i)
1430 {
1431 n = tree->a_nodes[i];
1432 dt = FABS(tree->times->nd_t[n->num]-tree->times->nd_t[n->anc->num]);
1433 n_jps = tree->times->n_jps[n->num];
1434 dens += Dpois(n_jps,lexp*dt,YES);
1435 }
1436
1437 tree->times->c_lnL_jps = dens;
1438
1439 return dens;
1440 }
1441
1442 //////////////////////////////////////////////////////////////
1443 //////////////////////////////////////////////////////////////
1444
1445
RATES_Posterior_Rates(t_tree * tree)1446 void RATES_Posterior_Rates(t_tree *tree)
1447 {
1448 int node_num;
1449 node_num = Rand_Int(0,2*tree->n_otu-3);
1450 RATES_Posterior_One_Rate(tree->a_nodes[node_num]->anc,tree->a_nodes[node_num],NO,tree);
1451 }
1452
1453 //////////////////////////////////////////////////////////////
1454 //////////////////////////////////////////////////////////////
1455
1456
RATES_Posterior_Times(t_tree * tree)1457 void RATES_Posterior_Times(t_tree *tree)
1458 {
1459 int node_num;
1460 node_num = Rand_Int(tree->n_otu,2*tree->n_otu-3);
1461 RATES_Posterior_One_Time(tree->a_nodes[node_num]->anc,tree->a_nodes[node_num],NO,tree);
1462 }
1463
1464 //////////////////////////////////////////////////////////////
1465 //////////////////////////////////////////////////////////////
1466
1467
RATES_Posterior_One_Rate(t_node * a,t_node * d,int traversal,t_tree * tree)1468 void RATES_Posterior_One_Rate(t_node *a, t_node *d, int traversal, t_tree *tree)
1469 {
1470 phydbl like_mean, like_var;
1471 phydbl prior_mean, prior_var;
1472 phydbl post_mean, post_var, post_sd;
1473 phydbl dt,rd,cr,cel,cvl;
1474 int dim;
1475 phydbl l_opp; /* length of the branch connected to the root, opposite to the one connected to d */
1476 t_edge *b;
1477 int i;
1478 t_node *v2,*v3;
1479 phydbl T0,T1,T2,T3;
1480 phydbl U0,U1,U2,U3;
1481 phydbl V1,V2,V3;
1482 phydbl r_min,r_max;
1483 int err;
1484 phydbl ratio,u;
1485 phydbl cvr, cer;
1486 phydbl nf;
1487 phydbl new_lnL_data, cur_lnL_data, new_lnL_rate, cur_lnL_rate;
1488 phydbl sd1,sd2,sd3;
1489 phydbl inflate_var;
1490
1491 if(d == tree->n_root) return;
1492
1493 dim = 2*tree->n_otu-3;
1494 err = NO;
1495
1496 inflate_var = tree->rates->inflate_var;
1497
1498 b = NULL;
1499 if(a == tree->n_root) b = tree->e_root;
1500 else for(i=0;i<3;i++) if(d->v[i] == a) { b = d->b[i]; break; }
1501
1502 v2 = v3 = NULL;
1503 if(!d->tax)
1504 {
1505 for(i=0;i<3;i++)
1506 if((d->v[i] != a) && (d->b[i] != tree->e_root))
1507 {
1508 if(!v2) { v2 = d->v[i]; }
1509 else { v3 = d->v[i]; }
1510 }
1511 }
1512
1513 T3 = T2 = 0.0;
1514 T0 = tree->times->nd_t[a->num];
1515 T1 = tree->times->nd_t[d->num];
1516 U0 = tree->rates->br_r[a->num];
1517 U1 = tree->rates->br_r[d->num];
1518 U3 = U2 = -1.0;
1519
1520 if(!d->tax)
1521 {
1522 T2 = tree->times->nd_t[v2->num];
1523 T3 = tree->times->nd_t[v3->num];
1524 U2 = tree->rates->br_r[v2->num];
1525 U3 = tree->rates->br_r[v3->num];
1526 }
1527
1528
1529 V1 = tree->rates->nu * FABS(T1 - T0);
1530 V2 = tree->rates->nu * FABS(T2 - T1);
1531 V3 = tree->rates->nu * FABS(T3 - T1);
1532
1533 dt = T1 - T0;
1534 rd = U1;
1535 cr = tree->rates->clock_r;
1536 r_min = tree->rates->min_rate;
1537 r_max = tree->rates->max_rate;
1538 nf = tree->rates->norm_fact;
1539 sd1 = SQRT(V1);
1540 sd2 = SQRT(V2);
1541 sd3 = SQRT(V3);
1542
1543
1544 /* Likelihood */
1545 cel=0.0;
1546 for(i=0;i<dim;i++)
1547 if(i != b->num)
1548 cel += tree->rates->reg_coeff[b->num*dim+i] * (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]);
1549 cel += tree->rates->mean_l[b->num];
1550 cvl = tree->rates->cond_var[b->num];
1551
1552 l_opp = 0.0;
1553 if(a == tree->n_root)
1554 {
1555 if(d == a->v[0]) l_opp = tree->rates->cur_l[a->v[1]->num];
1556 else l_opp = tree->rates->cur_l[a->v[0]->num];
1557 cel -= l_opp;
1558 }
1559
1560
1561 if(isnan(cvl) || isnan(cel) || cvl < .0)
1562 {
1563 for(i=0;i<dim;i++) if(i != b->num) printf("\n. reg: %f %f %f nu=%f clock=%f",
1564 tree->rates->reg_coeff[b->num*dim+i],
1565 tree->rates->u_cur_l[i],
1566 tree->rates->mean_l[i],
1567 tree->rates->nu,tree->rates->clock_r);
1568 PhyML_Fprintf(stderr,"\n. cel=%f cvl=%f\n",cel,cvl);
1569 PhyML_Fprintf(stderr,"\n. Warning: invalid expected and/or std. dev. values. Skipping this step.\n");
1570 Exit("\n");
1571 }
1572
1573 /* Model rates */
1574 /* if(tree->mod->log_l == YES) cel = exp(cel); */
1575 /* like_mean = cel / (dt*cr*nf); */
1576 /* like_var = cvl / POW(dt*cr*nf,2); */
1577
1578 like_mean = cel / (dt*cr*nf);
1579 like_var = cvl / POW(dt*cr*nf,2);
1580
1581 /* Prior */
1582 if(!d->tax)
1583 {
1584 cvr = 1./(1./V1 + 1./V2 + 1./V3);
1585 cer = cvr*(U0/V1 + U2/V2 + U3/V3);
1586 }
1587 else
1588 {
1589 cvr = V1;
1590 cer = U0;
1591 }
1592
1593 if(cvr < 0.0)
1594 {
1595 PhyML_Fprintf(stderr,"\n. cvr=%f d->tax=%d V1=%f v2=%f V3=%f",cvr,d->tax,V1,V2,V3);
1596 PhyML_Fprintf(stderr,"\n. T0=%f T1=%f T2=%f T3=%f",T0,T1,T2,T3);
1597 Exit("\n");
1598 }
1599
1600 prior_mean = cer;
1601 prior_var = cvr;
1602
1603 /* Posterior */
1604 post_mean = (prior_mean/prior_var + like_mean/like_var)/(1./prior_var + 1./like_var);
1605 post_var = 1./(1./prior_var + 1./like_var);
1606
1607
1608 /* Sample according to priors */
1609 if(tree->eval_alnL == NO)
1610 {
1611 post_mean = prior_mean;
1612 post_var = prior_var;
1613 }
1614
1615 post_sd = SQRT(post_var);
1616
1617 rd = Rnorm_Trunc(post_mean,inflate_var*post_sd,r_min,r_max,&err);
1618
1619 if(err || isnan(rd))
1620 {
1621 PhyML_Fprintf(stderr,"\n");
1622 PhyML_Fprintf(stderr,"\n. run: %d err=%d d->tax=%d",tree->mcmc->run,err,d->tax);
1623 PhyML_Fprintf(stderr,"\n. rd=%f cvr=%G dt=%G cr=%G",rd,cvr,dt,cr);
1624 PhyML_Fprintf(stderr,"\n. prior_mean = %G prior_var = %G",prior_mean,prior_var);
1625 PhyML_Fprintf(stderr,"\n. like_mean = %G like_var = %G",like_mean,like_var);
1626 PhyML_Fprintf(stderr,"\n. post_mean = %G post_var = %G",post_mean,post_var);
1627 PhyML_Fprintf(stderr,"\n. clock_r = %f",tree->rates->clock_r);
1628 PhyML_Fprintf(stderr,"\n. T0=%f T1=%f T2=%f T3=%f",T0,T1,T2,T3);
1629 PhyML_Fprintf(stderr,"\n. U0=%f U1=%f U2=%f U3=%f",U0,U1,U2,U3);
1630 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1631 Exit("\n");
1632 }
1633
1634 /* !!!!!!!!!!!!!! */
1635 /* u = Uni(); */
1636 /* rd = U1 * exp(1.*(u-0.5)); */
1637
1638 if(rd > r_min && rd < r_max)
1639 {
1640
1641 cur_lnL_data = tree->c_lnL;
1642 cur_lnL_rate = tree->rates->c_lnL_rates;
1643 new_lnL_data = tree->c_lnL;
1644 new_lnL_rate = tree->rates->c_lnL_rates;
1645
1646 tree->rates->br_r[d->num] = rd;
1647 RATES_Update_Norm_Fact(tree);
1648 RATES_Update_Cur_Bl(tree);
1649
1650 if(tree->eval_alnL) new_lnL_data = Lk(b,tree);
1651 /* new_lnL_rate = RATES_Lk_Rates(tree); */
1652 new_lnL_rate =
1653 cur_lnL_rate -
1654 (Log_Dnorm_Trunc(U1,U0,sd1,r_min,r_max,&err)) +
1655 (Log_Dnorm_Trunc(rd,U0,sd1,r_min,r_max,&err));
1656
1657 if(!d->tax)
1658 {
1659 new_lnL_rate -= (Log_Dnorm_Trunc(U2,U1,sd2,r_min,r_max,&err) + Log_Dnorm_Trunc(U3,U1,sd3,r_min,r_max,&err));
1660 new_lnL_rate += (Log_Dnorm_Trunc(U2,rd,sd2,r_min,r_max,&err) + Log_Dnorm_Trunc(U3,rd,sd3,r_min,r_max,&err));
1661 }
1662
1663 tree->rates->c_lnL_rates = new_lnL_rate;
1664
1665 /* printf("\n. %f %f sd1=%f U1=%f rd=%f ra=%f a=%d d=%d [%f] [%f %f]", */
1666 /* new_lnL_rate,RATES_Lk_Rates(tree),sd1,U1,rd,ra,a->num,d->num,tree->rates->br_r[tree->n_root->num], */
1667 /* Log_Dnorm_Trunc(U1,ra,sd1,r_min,r_max,&err), */
1668 /* Log_Dnorm_Trunc(rd,ra,sd1,r_min,r_max,&err)); */
1669
1670 ratio = 0.0;
1671 /* Proposal ratio */
1672 ratio += (Log_Dnorm_Trunc(U1,post_mean,inflate_var*post_sd,r_min,r_max,&err) - Log_Dnorm_Trunc(rd,post_mean,inflate_var*post_sd,r_min,r_max,&err));
1673 /* ratio += log(rd/U1); */
1674 /* Prior ratio */
1675 ratio += (new_lnL_rate - cur_lnL_rate);
1676 /* Likelihood ratio */
1677 ratio += (new_lnL_data - cur_lnL_data);
1678
1679
1680 ratio = exp(ratio);
1681
1682 /* printf("\n. R a=%3d T0=%6.1f T1=%6.1f T2=%6.1f T3=%6.1f ratio=%8f pm=%7f U1=%7.2f rd=%7.2f %f %f lr=%f %f ld=%f %f [%f]",a->num,T0,T1,T2,T3,ratio,post_mean,U1,rd, */
1683 /* Log_Dnorm_Trunc(U1,post_mean,post_sd,r_min,r_max,&err), */
1684 /* Log_Dnorm_Trunc(rd,post_mean,post_sd,r_min,r_max,&err), */
1685 /* new_lnL_rate,cur_lnL_rate, */
1686 /* new_lnL_data,cur_lnL_data, */
1687 /* ((Pnorm(r_max,U1,sd2)-Pnorm(r_min,U1,sd2)) * */
1688 /* (Pnorm(r_max,U1,sd3)-Pnorm(r_min,U1,sd3)))/ */
1689 /* ((Pnorm(r_max,rd,sd2)-Pnorm(r_min,rd,sd2)) * */
1690 /* (Pnorm(r_max,rd,sd3)-Pnorm(r_min,rd,sd3)))); */
1691
1692
1693 u = Uni();
1694
1695 if(u > MIN(1.,ratio))
1696 {
1697 tree->rates->br_r[d->num] = U1; /* reject */
1698 tree->rates->c_lnL_rates = cur_lnL_rate;
1699 tree->c_lnL = cur_lnL_data;
1700 RATES_Update_Cur_Bl(tree);
1701 Update_PMat_At_Given_Edge(b,tree);
1702 }
1703 else
1704 {
1705 tree->mcmc->acc_move[tree->mcmc->num_move_br_r+d->num]++;
1706 }
1707
1708 RATES_Update_Norm_Fact(tree);
1709 tree->mcmc->acc_rate[tree->mcmc->num_move_br_r+d->num] =
1710 (tree->mcmc->acc_move[tree->mcmc->num_move_br_r+d->num]+1.E-6)/
1711 (tree->mcmc->run_move[tree->mcmc->num_move_br_r+d->num]+1.E-6);
1712 }
1713
1714 tree->mcmc->run_move[tree->mcmc->num_move_br_r+d->num]++;
1715
1716
1717 if(traversal == YES)
1718 {
1719 if(d->tax == YES) return;
1720 else
1721 {
1722 for(i=0;i<3;i++)
1723 if(d->v[i] != a && d->b[i] != tree->e_root)
1724 {
1725 if(tree->io->lk_approx == EXACT) Update_Partial_Lk(tree,d->b[i],d);
1726 /* if(tree->io->lk_approx == EXACT) { tree->both_sides = YES; Lk(tree); } */
1727 RATES_Posterior_One_Rate(d,d->v[i],YES,tree);
1728 }
1729 }
1730
1731 if(tree->io->lk_approx == EXACT) Update_Partial_Lk(tree,b,d);
1732 /* if(tree->io->lk_approx == EXACT) { tree->both_sides = YES; Lk(tree); } */
1733 }
1734 }
1735
1736 //////////////////////////////////////////////////////////////
1737 //////////////////////////////////////////////////////////////
1738
1739
RATES_Posterior_One_Time(t_node * a,t_node * d,int traversal,t_tree * tree)1740 void RATES_Posterior_One_Time(t_node *a, t_node *d, int traversal, t_tree *tree)
1741 {
1742 /*
1743 T0 (a)
1744 |
1745 | l1,U1,b1
1746 |
1747 T1 (d)
1748 / \
1749 l2,u2,b2 / \ L3,u3,b3
1750 / \
1751 t2 t3
1752 (v2) (v3)
1753
1754 */
1755
1756 phydbl l1,l2,l3;
1757 phydbl new_l1;
1758 phydbl El1,El2,El3;
1759 phydbl u0,r1,r2,r3;
1760 phydbl t0,t1,t2,t3;
1761 phydbl t_min, t_max;
1762 /* phydbl t_max_12,t_max_13; */
1763 phydbl bl_min, bl_max;
1764 phydbl t1_new;
1765 phydbl EX,EY;
1766 phydbl VX,VY;
1767 phydbl cr;
1768 int i,j;
1769 phydbl *mu, *cov;
1770 phydbl *cond_mu, *cond_cov;
1771 short int *is_1;
1772 phydbl sig11, sig1X, sig1Y, sigXX, sigXY, sigYY;
1773 phydbl cov11,cov12,cov13,cov22,cov23,cov33;
1774 int dim;
1775 t_edge *b1, *b2, *b3;
1776 t_node *v2,*v3;
1777 phydbl *l2XY;
1778 phydbl l_opp;
1779 t_edge *buff_b;
1780 t_node *buff_n;
1781 int err;
1782 int num_1, num_2, num_3;
1783 phydbl nf;
1784 phydbl u, ratio;
1785 phydbl new_lnL_data, cur_lnL_data, new_lnL_rate, cur_lnL_rate;
1786 int num_move;
1787 phydbl inflate_var;
1788
1789 dim = 2*tree->n_otu-3;
1790 num_move = tree->mcmc->num_move_times;
1791 inflate_var = tree->rates->inflate_var;
1792
1793 if(d->tax) return;
1794
1795 if(FABS(tree->times->t_prior_min[d->num] - tree->times->t_prior_max[d->num]) < 1.E-10) return;
1796
1797 l2XY = tree->rates->_2n_vect2;
1798 mu = tree->rates->_2n_vect3;
1799 cov = tree->rates->_2n2n_vect1;
1800 cond_mu = tree->rates->_2n_vect1;
1801 cond_cov = tree->rates->_2n2n_vect2;
1802 is_1 = tree->rates->_2n_vect5;
1803 err = NO;
1804
1805 b1 = NULL;
1806 if(a == tree->n_root) b1 = tree->e_root;
1807 else for(i=0;i<3;i++) if(d->v[i] == a) { b1 = d->b[i]; break; }
1808
1809 b2 = b3 = NULL;
1810 v2 = v3 = NULL;
1811 for(i=0;i<3;i++)
1812 if((d->v[i] != a) && (d->b[i] != tree->e_root))
1813 {
1814 if(!v2) { v2 = d->v[i]; b2 = d->b[i]; }
1815 else { v3 = d->v[i]; b3 = d->b[i]; }
1816 }
1817
1818 t2 = tree->times->nd_t[v2->num];
1819 t3 = tree->times->nd_t[v3->num];
1820
1821 buff_n = NULL;
1822 buff_b = NULL;
1823 if(t3 > t2)
1824 {
1825 buff_n = v2;
1826 v2 = v3;
1827 v3 = buff_n;
1828
1829 buff_b = b2;
1830 b2 = b3;
1831 b3 = buff_b;
1832 }
1833
1834 t0 = tree->times->nd_t[a->num];
1835 t1 = tree->times->nd_t[d->num];
1836 t2 = tree->times->nd_t[v2->num];
1837 t3 = tree->times->nd_t[v3->num];
1838 u0 = tree->rates->br_r[a->num];
1839 r1 = tree->rates->br_r[d->num];
1840 r2 = tree->rates->br_r[v2->num];
1841 r3 = tree->rates->br_r[v3->num];
1842 l1 = tree->rates->cur_l[d->num];
1843 l2 = tree->rates->cur_l[v2->num];
1844 l3 = tree->rates->cur_l[v3->num];
1845 cr = tree->rates->clock_r;
1846 nf = tree->rates->norm_fact;
1847
1848 for(i=0;i<dim;i++) is_1[i] = 0;
1849 is_1[b1->num] = 1;
1850 is_1[b2->num] = 1;
1851 is_1[b3->num] = 1;
1852
1853 /* for(i=0;i<dim;i++) cond_mu[i] = 0.0; */
1854 /* For(i,dim*dim) cond_cov[i] = 0.0; */
1855 /* Normal_Conditional(tree->rates->mean_l,tree->rates->cov_l,tree->rates->u_cur_l,dim,is_1,3,cond_mu,cond_cov); */
1856
1857 /* El1 = cond_mu[b1->num]; */
1858 /* El2 = cond_mu[b2->num]; */
1859 /* El3 = cond_mu[b3->num]; */
1860
1861 /* cov11 = cond_cov[b1->num*dim+b1->num]; */
1862 /* cov12 = cond_cov[b1->num*dim+b2->num]; */
1863 /* cov13 = cond_cov[b1->num*dim+b3->num]; */
1864 /* cov23 = cond_cov[b2->num*dim+b3->num]; */
1865 /* cov22 = cond_cov[b2->num*dim+b2->num]; */
1866 /* cov33 = cond_cov[b3->num*dim+b3->num]; */
1867
1868
1869 /* El1 = tree->rates->u_ml_l[b1->num]; */
1870 /* El2 = tree->rates->u_ml_l[b2->num]; */
1871 /* El3 = tree->rates->u_ml_l[b3->num]; */
1872
1873 /* cov11 = tree->rates->cov[b1->num*dim+b1->num]; */
1874 /* cov12 = tree->rates->cov[b1->num*dim+b2->num]; */
1875 /* cov13 = tree->rates->cov[b1->num*dim+b3->num]; */
1876 /* cov23 = tree->rates->cov[b2->num*dim+b3->num]; */
1877 /* cov22 = tree->rates->cov[b2->num*dim+b2->num]; */
1878 /* cov33 = tree->rates->cov[b3->num*dim+b3->num]; */
1879
1880 /* PhyML_Fprintf(stderr,"\n- El1=%10f El2=%10f El3=%10f",El1,El2,El3); */
1881 /* PhyML_Fprintf(stderr,"\n- cov11=%10f cov12=%10f cov13=%10f cov23=%10f cov22=%10f cov33=%10f", cov11,cov12,cov13,cov23,cov22,cov33); */
1882
1883 num_1 = num_2 = num_3 = -1;
1884 if(b1->num < b2->num && b2->num < b3->num) { num_1 = 0; num_2 = 1; num_3 = 2; }
1885 if(b1->num < b3->num && b3->num < b2->num) { num_1 = 0; num_2 = 2; num_3 = 1; }
1886 if(b2->num < b1->num && b1->num < b3->num) { num_1 = 1; num_2 = 0; num_3 = 2; }
1887 if(b2->num < b3->num && b3->num < b1->num) { num_1 = 2; num_2 = 0; num_3 = 1; }
1888 if(b3->num < b2->num && b2->num < b1->num) { num_1 = 2; num_2 = 1; num_3 = 0; }
1889 if(b3->num < b1->num && b1->num < b2->num) { num_1 = 1; num_2 = 2; num_3 = 0; }
1890
1891 cov11 = tree->rates->trip_cond_cov[d->num * 9 + num_1 * 3 + num_1];
1892 cov12 = tree->rates->trip_cond_cov[d->num * 9 + num_1 * 3 + num_2];
1893 cov13 = tree->rates->trip_cond_cov[d->num * 9 + num_1 * 3 + num_3];
1894 cov23 = tree->rates->trip_cond_cov[d->num * 9 + num_2 * 3 + num_3];
1895 cov22 = tree->rates->trip_cond_cov[d->num * 9 + num_2 * 3 + num_2];
1896 cov33 = tree->rates->trip_cond_cov[d->num * 9 + num_3 * 3 + num_3];
1897
1898 El1=0.0;
1899 for(i=0;i<dim;i++)
1900 if(i != b1->num && i != b2->num && i != b3->num)
1901 El1 += tree->rates->trip_reg_coeff[d->num * (6*tree->n_otu-9) + num_1 * dim +i] * (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]);
1902 El1 += tree->rates->mean_l[b1->num];
1903
1904 El2=0.0;
1905 for(i=0;i<dim;i++)
1906 if(i != b1->num && i != b2->num && i != b3->num)
1907 El2 += tree->rates->trip_reg_coeff[d->num * (6*tree->n_otu-9) + num_2 * dim +i] * (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]);
1908 El2 += tree->rates->mean_l[b2->num];
1909
1910 El3=0.0;
1911 for(i=0;i<dim;i++)
1912 if(i != b1->num && i != b2->num && i != b3->num)
1913 El3 += tree->rates->trip_reg_coeff[d->num * (6*tree->n_otu-9) + num_3 * dim +i] * (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]);
1914 El3 += tree->rates->mean_l[b3->num];
1915
1916
1917 /* PhyML_Fprintf(stderr,"\n+ El1=%10f El2=%10f El3=%10f",El1,El2,El3); */
1918 /* PhyML_Fprintf(stderr,"\n+ cov11=%10f cov12=%10f cov13=%10f cov23=%10f cov22=%10f cov33=%10f", cov11,cov12,cov13,cov23,cov22,cov33); */
1919
1920
1921 t1_new = +1;
1922
1923 t_min = MAX(t0,tree->times->t_prior_min[d->num]);
1924 t_max = MIN(MIN(t2,t3),tree->times->t_prior_max[d->num]);
1925
1926 t_min += tree->rates->min_dt;
1927 t_max -= tree->rates->min_dt;
1928
1929 if(t_min > t_max)
1930 {
1931 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1932 Exit("\n");
1933 }
1934
1935 bl_min = bl_max = -1.0;
1936
1937 l_opp = 0.0;
1938 if(a == tree->n_root)
1939 {
1940 l_opp = (d == a->v[0])?(tree->rates->cur_l[a->v[1]->num]):(tree->rates->cur_l[a->v[0]->num]);
1941 El1 -= l_opp;
1942 }
1943
1944 EX = El1/r1 + El2/r2;
1945 EY = El1/r1 + El3/r3;
1946
1947 VX = cov11/(r1*r1) + cov22/(r2*r2) + 2.*cov12/(r1*r2);
1948 VY = cov11/(r1*r1) + cov33/(r3*r3) + 2.*cov13/(r1*r3);
1949
1950 mu[0] = El1;
1951 mu[1] = EX;
1952 mu[2] = EY;
1953
1954 sig11 = cov11;
1955 sig1X = cov11/r1 + cov12/r2;
1956 sig1Y = cov11/r1 + cov13/r3;
1957 sigXX = VX;
1958 sigYY = VY;
1959 sigXY = cov11/(r1*r1) + cov13/(r1*r3) + cov12/(r1*r2) + cov23/(r2*r3);
1960
1961 cov[0*3+0] = sig11; cov[0*3+1] = sig1X; cov[0*3+2] = sig1Y;
1962 cov[1*3+0] = sig1X; cov[1*3+1] = sigXX; cov[1*3+2] = sigXY;
1963 cov[2*3+0] = sig1Y; cov[2*3+1] = sigXY; cov[2*3+2] = sigYY;
1964
1965 l2XY[0] = 0.0; /* does not matter */
1966 l2XY[1] = (t2-t0)*cr*nf; /* constraint 1 */
1967 l2XY[2] = (t3-t0)*cr*nf; /* constraint 2 */
1968
1969 is_1[0] = is_1[1] = is_1[2] = 0;
1970 is_1[0] = 1;
1971
1972 Normal_Conditional(mu,cov,l2XY,3,is_1,1,cond_mu,cond_cov);
1973
1974
1975 if(cond_cov[0*3+0] < 0.0)
1976 {
1977 PhyML_Printf("\n. a: %d d: %d",a->num,d->num);
1978 PhyML_Printf("\n. Conditional mean=%G var=%G",cond_mu[0],cond_cov[0*3+0]);
1979 PhyML_Printf("\n. t0=%G t1=%f t2=%f t3=%f l1=%G l2=%G l3=%G",t0,t1,t2,t3,l1,l2,l3);
1980 PhyML_Printf("\n. El1=%G El2=%G El3=%G Nu=%G r1=%G r2=%G r3=%G cr=%G",El1,El2,El3,tree->rates->nu,r1,r2,r3,cr);
1981 PhyML_Printf("\n. COV11=%f COV12=%f COV13=%f COV22=%f COV23=%f COV33=%f",cov11,cov12,cov13,cov22,cov23,cov33);
1982 PhyML_Printf("\n. constraint1: %f constraints2: %f",l2XY[1],l2XY[2]);
1983 PhyML_Printf("\n");
1984 for(i=0;i<3;i++)
1985 {
1986 PhyML_Printf(". mu%d=%12lf\t",i,mu[i]);
1987 for(j=0;j<3;j++)
1988 {
1989 PhyML_Printf("%12lf ",cov[i*3+j]);
1990 }
1991 PhyML_Printf("\n");
1992 }
1993 cond_cov[0*3+0] = 1.E-10;
1994 }
1995
1996
1997 bl_min = (t_min - t0) * r1 * cr * nf;
1998 bl_max = (t_max - t0) * r1 * cr * nf;
1999
2000 new_l1 = Rnorm_Trunc(cond_mu[0],inflate_var*SQRT(cond_cov[0*3+0]),bl_min,bl_max,&err);
2001 /* new_l1 = Rnorm(cond_mu[0],SQRT(cond_cov[0*3+0])); */
2002
2003
2004 if(new_l1 < bl_min) new_l1 = l1;
2005 if(new_l1 > bl_max) new_l1 = l1;
2006
2007 t1_new = new_l1/(r1*cr*nf) + t0;
2008
2009
2010 if(err)
2011 {
2012 PhyML_Printf("\n");
2013 PhyML_Printf("\n. Root ? %s",(tree->n_root==a)?("yes"):("no"));
2014 PhyML_Printf("\n. %s %d %d",__FILE__,__LINE__,tree->mcmc->run);
2015 PhyML_Printf("\n. t0=%f t1=%f t2=%f t3=%f",t0,t1,t2,t3);
2016 PhyML_Printf("\n. t_min=%f t_max=%f",t_min,t_max);
2017 PhyML_Printf("\n. bl_min=%f bl_max=%f",bl_min,bl_max);
2018 PhyML_Printf("\n. cond_mu[0]=%f cond_cov[0]=%f",cond_mu[0],SQRT(cond_cov[0*3+0]));
2019 PhyML_Printf("\n. El1=%f El2=%f El3=%f",El1,El2,El3);
2020 PhyML_Printf("\n. l1=%f l2=%f l3=%f",l1,l2,l3);
2021 PhyML_Printf("\n. u0=%f r1=%f r2=%f r3=%f",u0,r1,r2,r3);
2022 PhyML_Printf("\n. COV11=%f COV22=%f",cov11,cov22);
2023 PhyML_Printf("\n. Clock rate = %f",tree->rates->clock_r);
2024 PhyML_Printf("\n. Setting t1_new to %f",t1);
2025 t1_new = t1;
2026 /* Exit("\n"); */
2027 }
2028 /* } */
2029 /* else */
2030 /* { */
2031 /* bl_min = (t2 - t_max) * r2; */
2032 /* bl_max = (t2 - t_min) * r2; */
2033
2034 /* l2 = Rnorm_Trunc(cond_mu[0],SQRT(cond_cov[0*3+0]),bl_min,bl_max,&err); */
2035 /* t1_new = -l2/r2 + t2; */
2036
2037 /* if(err) */
2038 /* { */
2039 /* PhyML_Printf("\n"); */
2040 /* PhyML_Printf("\n. %s %d %d",__FILE__,__LINE__,tree->mcmc->run); */
2041 /* PhyML_Printf("\n. t0=%f t1=%f t2=%f t3=%f",t0,t1,t2,t3); */
2042 /* PhyML_Printf("\n. t_min=%f t_max=%f",t_min,t_max); */
2043 /* PhyML_Printf("\n. bl_min=%f bl_max=%f",bl_min,bl_max); */
2044 /* PhyML_Printf("\n. cond_mu[0]=%f cond_cov[0]=%f",cond_mu[0],SQRT(cond_cov[0*3+0])); */
2045 /* PhyML_Printf("\n. El1=%f El2=%f El3=%f",El1,El2,El3); */
2046 /* PhyML_Printf("\n. l1=%f l2=%f l3=%f",l1,l2,l3); */
2047 /* PhyML_Printf("\n. u0=%f r1=%f r2=%f r3=%f",u0,r1,r2,r3); */
2048 /* PhyML_Printf("\n. COV11=%f COV22=%f",cov11,cov22); */
2049 /* PhyML_Printf("\n. Clock rate = %f",tree->rates->clock_r); */
2050 /* PhyML_Printf("\n. Setting t1_new to %f",t1); */
2051 /* t1_new = t1; */
2052 /* /\* Exit("\n"); *\/ */
2053 /* } */
2054 /* } */
2055
2056 if(t1_new < t0)
2057 {
2058 t1_new = t0+1.E-4;
2059 PhyML_Printf("\n");
2060 PhyML_Printf("\n. a is root -> %s",(a == tree->n_root)?("YES"):("NO"));
2061 PhyML_Printf("\n. t0 = %f t1_new = %f t1 = %f",t0,t1_new,t1);
2062 PhyML_Printf("\n. t_min=%f t_max=%f",t_min,t_max);
2063 PhyML_Printf("\n. l1 = %f",l1);
2064 PhyML_Printf("\n. bl_min = %f bl_max = %f",bl_min,bl_max);
2065 PhyML_Printf("\n. (t1-t0)=%f (t2-t1)=%f",t1-t0,t2-t1);
2066 PhyML_Printf("\n. l1 = %f l2 = %f cov11=%f cov22=%f cov33=%f",l1,l2,cov11,cov22,cov33);
2067 PhyML_Printf("\n. clock=%G",tree->rates->clock_r);
2068 PhyML_Printf("\n. u0=%f r1=%f r2=%f r3=%f",u0,r1,r2,r3);
2069 PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
2070 /* Exit("\n"); */
2071 }
2072 if(t1_new > MIN(t2,t3))
2073 {
2074 PhyML_Printf("\n");
2075 PhyML_Printf("\n. a is root -> %s",(a == tree->n_root)?("YES"):("NO"));
2076 PhyML_Printf("\n. t0 = %f t1_new = %f t1 = %f t2 = %f t3 = %f MIN(t2,t3)=%f",t0,t1_new,t1,t2,t3,MIN(t2,t3));
2077 PhyML_Printf("\n. t_min=%f t_max=%f",t_min,t_max);
2078 PhyML_Printf("\n. l2 = %f",l2);
2079 PhyML_Printf("\n. bl_min = %f bl_max = %f",bl_min,bl_max);
2080 PhyML_Printf("\n. (t1-t0)=%f (t2-t1)=%f",t1-t0,t2-t1);
2081 PhyML_Printf("\n. l1 = %f l2 = %f cov11=%f cov22=%f cov33=%f",l1,l2,cov11,cov22,cov33);
2082 PhyML_Printf("\n. clock=%G",tree->rates->clock_r);
2083 PhyML_Printf("\n. u0=%f r1=%f r2=%f r3=%f",u0,r1,r2,r3);
2084 PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
2085 /* Exit("\n"); */
2086 }
2087
2088 if(isnan(t1_new))
2089 {
2090 PhyML_Printf("\n. run=%d",tree->mcmc->run);
2091 PhyML_Printf("\n. mean=%G var=%G",cond_mu[0],cond_cov[0*3+0]);
2092 PhyML_Printf("\n. t1=%f l1=%G r1=%G t0=%G",t1,l1,r1,t0);
2093 PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
2094 /* Exit("\n"); */
2095 }
2096
2097 if(tree->eval_alnL == YES)
2098 tree->times->nd_t[d->num] = t1_new;
2099 else
2100 {
2101 u = Uni();
2102 tree->times->nd_t[d->num] = u*(t_max-t_min) + t_min;
2103 }
2104
2105 RATES_Update_Norm_Fact(tree);
2106 RATES_Update_Cur_Bl(tree);
2107
2108 cur_lnL_data = tree->c_lnL;
2109 cur_lnL_rate = tree->rates->c_lnL_rates;
2110 new_lnL_data = tree->c_lnL;
2111 new_lnL_rate = tree->rates->c_lnL_rates;
2112
2113 /* !!!!!!!!!!!!!!!!! */
2114 /* tree->both_sides = NO; */
2115 /* new_lnL_data = Lk(tree); */
2116
2117 if(tree->io->lk_approx == EXACT)
2118 {
2119 Update_PMat_At_Given_Edge(b1,tree);
2120 Update_PMat_At_Given_Edge(b2,tree);
2121 Update_PMat_At_Given_Edge(b3,tree);
2122 Update_Partial_Lk(tree,b1,d);
2123 }
2124 new_lnL_data = Lk(b1,tree);
2125
2126 new_lnL_rate = RATES_Lk_Rates(tree);
2127
2128 ratio = 0.0;
2129
2130 /* Proposal ratio */
2131 if(tree->eval_alnL)
2132 ratio += (Log_Dnorm_Trunc(l1, cond_mu[0],inflate_var*SQRT(cond_cov[0*3+0]),bl_min,bl_max,&err) -
2133 Log_Dnorm_Trunc(new_l1,cond_mu[0],inflate_var*SQRT(cond_cov[0*3+0]),bl_min,bl_max,&err));
2134
2135 /* Prior ratio */
2136 ratio += (new_lnL_rate - cur_lnL_rate);
2137
2138 /* Likelihood ratio */
2139 ratio += (new_lnL_data - cur_lnL_data);
2140
2141 /* printf("\n* d:%d Ratio=%f l1=%f new_l1=%f mean=%f ml=%f sd=%f [%f %f]", */
2142 /* d->num, */
2143 /* exp(ratio), */
2144 /* l1,new_l1, */
2145 /* cond_mu[0], */
2146 /* tree->rates->mean_l[b1->num], */
2147 /* SQRT(cond_cov[0*3+0]), */
2148 /* Log_Dnorm(l1,cond_mu[0],SQRT(cond_cov[0*3+0]),&err),Log_Dnorm(new_l1,cond_mu[0],SQRT(cond_cov[0*3+0]),&err)); */
2149
2150 ratio = exp(ratio);
2151
2152
2153 u = Uni();
2154 if(u > MIN(1.,ratio))
2155 {
2156 tree->times->nd_t[d->num] = t1; /* reject */
2157 tree->rates->c_lnL_rates = cur_lnL_rate;
2158 tree->c_lnL = cur_lnL_data;
2159 RATES_Update_Cur_Bl(tree);
2160 if(tree->io->lk_approx == EXACT)
2161 {
2162 Update_PMat_At_Given_Edge(b1,tree);
2163 Update_PMat_At_Given_Edge(b2,tree);
2164 Update_PMat_At_Given_Edge(b3,tree);
2165 Update_Partial_Lk(tree,b1,d);
2166 }
2167 }
2168 else
2169 {
2170 tree->mcmc->acc_move[num_move]++;
2171 }
2172
2173 RATES_Update_Norm_Fact(tree);
2174
2175 tree->mcmc->run_move[num_move]++;
2176
2177 if(traversal == YES)
2178 {
2179 if(d->tax == YES) return;
2180 else
2181 {
2182 for(i=0;i<3;i++)
2183 if(d->v[i] != a && d->b[i] != tree->e_root)
2184 {
2185 if(tree->io->lk_approx == EXACT) Update_Partial_Lk(tree,d->b[i],d);
2186 RATES_Posterior_One_Time(d,d->v[i],YES,tree);
2187 }
2188 }
2189 if(tree->io->lk_approx == EXACT) Update_Partial_Lk(tree,b1,d);
2190 }
2191
2192 }
2193
2194 //////////////////////////////////////////////////////////////
2195 //////////////////////////////////////////////////////////////
2196
2197
RATES_Posterior_Time_Root(t_tree * tree)2198 void RATES_Posterior_Time_Root(t_tree *tree)
2199 {
2200 /*
2201 Root, T0
2202 / \
2203 l1 / \ l2
2204 / \
2205 v1,t1 v2,t2
2206
2207 */
2208
2209 phydbl t1,t2;
2210 phydbl u0,u1,u2;
2211 phydbl cel,cvl;
2212 int i,dim;
2213 t_edge *b;
2214 t_node *root;
2215 phydbl t0,t0_min, t0_max;
2216 phydbl new_t0;
2217 /* phydbl t_max_01, t_max_02; */
2218 int err;
2219 phydbl cr;
2220 phydbl bl_min, bl_max;
2221 phydbl new_l;
2222 phydbl new_lnL_data, cur_lnL_data, cur_lnL_rate;
2223 phydbl u,ratio;
2224
2225 dim = 2*tree->n_otu-3;
2226 b = tree->e_root;
2227 root = tree->n_root;
2228 t0 = tree->times->nd_t[root->num];
2229 t1 = tree->times->nd_t[root->v[2]->num];
2230 t2 = tree->times->nd_t[root->v[1]->num];
2231 u1 = tree->rates->br_r[root->v[2]->num];
2232 u2 = tree->rates->br_r[root->v[1]->num];
2233 u0 = 1.0;
2234 cr = tree->rates->clock_r;
2235 t0_min = -BIG;
2236 t0_max = MIN(t1,t2);
2237
2238 /* t0_max = MIN(t1 - (1./tree->rates->nu)*POW((u0-u1)/tree->rates->z_max,2), */
2239 /* t2 - (1./tree->rates->nu)*POW((u0-u2)/tree->rates->z_max,2)); */
2240
2241 /* if(u1 > u0) t_max_01 = t1 - (1./tree->rates->nu)*POW((u1-u0)/(PointNormal(tree->rates->p_max*Pnorm(.0,.0,1.))),2); */
2242 /* else t_max_01 = t1 - (1./tree->rates->nu)*POW((u1-u0)/(PointNormal(tree->rates->p_max*(1.-Pnorm(.0,.0,1.)))),2); */
2243
2244 /* if(u2 > u0) t_max_02 = t2 - (1./tree->rates->nu)*POW((u2-u0)/(PointNormal(tree->rates->p_max*Pnorm(.0,.0,1.))),2); */
2245 /* else t_max_02 = t2 - (1./tree->rates->nu)*POW((u2-u0)/(PointNormal(tree->rates->p_max*(1.-Pnorm(.0,.0,1.)))),2); */
2246
2247
2248 /* t_max_01 = RATES_Find_Max_Dt_Bisec(u1,u0,tree->times->t_prior_min[root->num],t1,tree->rates->nu,tree->rates->p_max,(u1 < u0)?YES:NO); */
2249 /* t_max_02 = RATES_Find_Max_Dt_Bisec(u2,u0,tree->times->t_prior_min[root->num],t2,tree->rates->nu,tree->rates->p_max,(u2 < u0)?YES:NO); */
2250 /* t0_max = MIN(t_max_01,t_max_02); */
2251
2252 /* RATES_Min_Max_Interval(u0,u1,u2,r3,t0,t2,t3,&t_min,&t_max,tree->rates->nu,tree->rates->p_max,tree); */
2253
2254 t0_min = MAX(t0_min,tree->times->t_prior_min[root->num]);
2255 t0_max = MIN(t0_max,tree->times->t_prior_max[root->num]);
2256
2257 t0_max -= tree->rates->min_dt;
2258
2259 if(t0_min > t0_max) return;
2260
2261 tree->times->t_prior[root->num] = Uni()*(t0_max - t0_min) + t0_min;
2262
2263 u0 *= cr;
2264 u1 *= cr;
2265 u2 *= cr;
2266
2267 if(FABS(tree->times->t_prior_min[root->num] - tree->times->t_prior_max[root->num]) < 1.E-10) return;
2268
2269 cel=0.0;
2270 for(i=0;i<dim;i++) if(i != b->num) cel += tree->rates->reg_coeff[b->num*dim+i] * (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]);
2271 cel += tree->rates->mean_l[b->num];
2272
2273 cvl = tree->rates->cond_var[b->num];
2274
2275 bl_min = u1 * (t1 - t0_max) + u2 * (t2 - t0_max);
2276 bl_max = u1 * (t1 - t0_min) + u2 * (t2 - t0_min);
2277
2278 if(bl_min > bl_max) return;
2279
2280 new_l = Rnorm_Trunc(cel,SQRT(cvl),bl_min,bl_max,&err);
2281
2282 new_t0 = (u1*t1 + u2*t2 - new_l)/(u1+u2);
2283
2284 if(t0 > t1 || t0 > t2)
2285 {
2286 PhyML_Printf("\n");
2287 PhyML_Printf("\n. Run = %d",tree->mcmc->run);
2288 PhyML_Printf("\n. t0=%f t1=%f t2=%f",t0,t1,t2);
2289 PhyML_Printf("\n. t0_min=%f t0_max=%f",t0_min,t0_max);
2290 PhyML_Printf("\n. new_l=%f cel=%f",new_l,cel);
2291 PhyML_Printf("\n. u0=%f u1=%f u2=%f",u0/cr,u1/cr,u2/cr);
2292 PhyML_Printf("\n. Nu = %f Clock = %f",tree->rates->nu,tree->rates->clock_r);
2293 PhyML_Printf("\n. Setting t0 to %f",tree->times->nd_t[root->num]);
2294 return;
2295 }
2296
2297 if(t0 < t0_min || t0 > t0_max)
2298 {
2299 PhyML_Printf("\n");
2300 PhyML_Printf("\n. Run = %d",tree->mcmc->run);
2301 PhyML_Printf("\n. t0=%f t1=%f t2=%f",t0,t1,t2);
2302 PhyML_Printf("\n. t0_min=%f t0_max=%f",t0_min,t0_max);
2303 PhyML_Printf("\n. u0=%f u1=%f u2=%f",u0/cr,u1/cr,u2/cr);
2304 PhyML_Printf("\n. Nu = %f Clock = %f",tree->rates->nu,tree->rates->clock_r);
2305 PhyML_Printf("\n. Setting t0 to %f",tree->times->nd_t[root->num]);
2306 t0 = tree->times->nd_t[root->num];
2307 }
2308
2309 /* Sample according to prior */
2310 /* tree->times->nd_t[root->num] = tree->times->t_prior[root->num]; */
2311
2312 /* Sample according to posterior */
2313 tree->times->nd_t[root->num] = new_t0;
2314
2315 RATES_Update_Norm_Fact(tree);
2316 RATES_Update_Cur_Bl(tree);
2317
2318 cur_lnL_data = tree->c_lnL;
2319 cur_lnL_rate = tree->rates->c_lnL_rates;
2320 new_lnL_data = tree->c_lnL;
2321
2322 new_lnL_data = Lk(NULL,tree);
2323
2324 ratio = 0.0;
2325 /* Prior ratio */
2326 ratio += .0;
2327 /* Likelihood ratio */
2328 ratio += new_lnL_data - cur_lnL_data;
2329
2330 ratio = exp(ratio);
2331 u = Uni();
2332 if(u > MIN(1.,ratio))
2333 {
2334 tree->times->nd_t[root->num] = t0; /* reject */
2335 tree->rates->c_lnL_rates = cur_lnL_rate;
2336 tree->c_lnL = cur_lnL_data;
2337 }
2338 else
2339 {
2340 tree->mcmc->acc_move[tree->mcmc->num_move_times]++;
2341 }
2342
2343 RATES_Update_Norm_Fact(tree);
2344 RATES_Update_Cur_Bl(tree);
2345
2346 tree->mcmc->run_move[tree->mcmc->num_move_times]++;
2347 tree->mcmc->acc_rate[tree->mcmc->num_move_times] =
2348 (tree->mcmc->acc_move[tree->mcmc->num_move_times]+1.E-6)/
2349 (tree->mcmc->run_move[tree->mcmc->num_move_times]+1.E-6);
2350
2351 }
2352
2353 //////////////////////////////////////////////////////////////
2354 //////////////////////////////////////////////////////////////
2355
RATES_Update_Cur_Bl(t_tree * tree)2356 void RATES_Update_Cur_Bl(t_tree *tree)
2357 {
2358
2359 RATES_Update_Cur_Bl_Pre(tree->n_root,tree->n_root->v[1],tree->n_root->b[1],tree);
2360 RATES_Update_Cur_Bl_Pre(tree->n_root,tree->n_root->v[2],tree->n_root->b[2],tree);
2361
2362 if(tree->mod && tree->mod->log_l == YES)
2363 {
2364 tree->e_root->l->v =
2365 exp(tree->n_root->b[1]->l->v) +
2366 exp(tree->n_root->b[2]->l->v) ;
2367 tree->e_root->l->v = log(tree->e_root->l->v);
2368 }
2369 else
2370 {
2371 tree->e_root->l->v = tree->n_root->b[1]->l->v + tree->n_root->b[2]->l->v;
2372 }
2373
2374 tree->rates->u_cur_l[tree->e_root->num] = tree->e_root->l->v;
2375 tree->n_root_pos = tree->n_root->b[2]->l->v / tree->e_root->l->v;
2376
2377 if(tree->rates->model == GUINDON)
2378 {
2379 phydbl t0,t1,t2;
2380 t_node *n0, *n1;
2381
2382 n0 = tree->n_root->v[2];
2383 n1 = tree->n_root->v[1];
2384 t1 = tree->times->nd_t[tree->n_root->v[2]->num];
2385 t2 = tree->times->nd_t[tree->n_root->v[1]->num];
2386 t0 = tree->times->nd_t[tree->n_root->num];
2387
2388 tree->e_root->l->v =
2389 (t1-t0)/(t1+t2-2.*t0)*tree->rates->cur_gamma_prior_mean[n0->num] +
2390 (t2-t0)/(t1+t2-2.*t0)*tree->rates->cur_gamma_prior_mean[n1->num];
2391
2392 tree->e_root->l_var->v =
2393 pow((t1-t0)/(t1+t2-2.*t0),2)*tree->rates->cur_gamma_prior_var[n0->num] +
2394 pow((t2-t0)/(t1+t2-2.*t0),2)*tree->rates->cur_gamma_prior_var[n1->num];
2395 }
2396
2397 if(tree->is_mixt_tree == YES) MIXT_RATES_Update_Cur_Bl(tree);
2398 }
2399
2400 //////////////////////////////////////////////////////////////
2401 //////////////////////////////////////////////////////////////
2402
RATES_Update_Cur_Bl_Pre(t_node * a,t_node * d,t_edge * b,t_tree * tree)2403 void RATES_Update_Cur_Bl_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
2404 {
2405 phydbl dt,rr,cr,ra,rd,ta,td,nu;
2406
2407 assert(a);
2408 assert(d);
2409
2410 ra = rd = -1.;
2411 tree->rates->br_do_updt[d->num] = YES;
2412
2413 if(tree->rates->br_do_updt[d->num] == YES)
2414 {
2415 tree->rates->br_do_updt[d->num] = NO;
2416
2417 if(tree->rates->model == LOGNORMAL ||
2418 tree->rates->model == THORNE ||
2419 tree->rates->model == STRICTCLOCK)
2420 {
2421 rd = tree->rates->br_r[d->num];
2422 ra = tree->rates->br_r[a->num];
2423 }
2424 else if(tree->rates->model == GUINDON)
2425 {
2426 rd = tree->rates->nd_r[d->num];
2427 ra = tree->rates->nd_r[a->num];
2428 }
2429 else assert(FALSE);
2430
2431 dt = fabs(tree->times->nd_t[d->num] - tree->times->nd_t[a->num]);
2432 cr = tree->rates->clock_r;
2433 td = tree->times->nd_t[d->num];
2434 ta = tree->times->nd_t[a->num];
2435 nu = tree->rates->nu;
2436 rr = -1.0;
2437
2438 if(tree->rates->model == LOGNORMAL)
2439 {
2440 rr = rd;
2441 tree->rates->cur_l[d->num] = dt*rr*cr;
2442 /* PhyML_Printf("\n. a: %d d: %d rr: %f dt: %f cr: %f tree: %p",a->num,d->num,rr,dt,cr,tree); */
2443 }
2444
2445 if(tree->rates->model == THORNE)
2446 {
2447 rr = (ra+rd)/2.;
2448 tree->rates->cur_l[d->num] = dt*rr*cr;
2449 }
2450
2451 if(tree->rates->model == GUINDON)
2452 {
2453 phydbl m,v;
2454
2455 Integrated_Geometric_Brownian_Bridge_Moments(dt,ra,rd,nu,&m,&v);
2456
2457 if(isnan(m) || isnan(v) || m < 0.0 || v < 0.0)
2458 {
2459 PhyML_Fprintf(stderr,"\n. dt: %G ra: %G rd: %G nu: %G m: %G v: %G a is root ? %d d is root ? %d",
2460 dt,ra,rd,nu,m,v,
2461 a==tree->n_root,
2462 d==tree->n_root);
2463 }
2464
2465 m *= cr*dt; // the actual rate average is m * cr. We multiply by dt in order to derive the value for the branch length
2466 v *= (cr*cr)*(dt*dt);
2467
2468 tree->rates->cur_gamma_prior_mean[d->num] = m;
2469 tree->rates->cur_gamma_prior_var[d->num] = v;
2470
2471 tree->rates->cur_l[d->num] = tree->rates->cur_gamma_prior_mean[d->num]; // Required for having proper branch lengths in Write_Tree function
2472 }
2473
2474 if(tree->rates->model == STRICTCLOCK)
2475 {
2476 tree->rates->cur_l[d->num] = dt*cr;
2477 }
2478
2479 /* printf("\n. td: %12f ta: %12f dt: %12f cr: %12f ra: %12f rd: %12f l: %12f", */
2480 /* tree->times->nd_t[d->num], */
2481 /* tree->times->nd_t[a->num], */
2482 /* dt,cr,ra,rd,tree->rates->cur_l[d->num]); */
2483
2484
2485 if(tree->mod && tree->mod->log_l == YES) tree->rates->cur_l[d->num] = log(tree->rates->cur_l[d->num]);
2486
2487 if(b)
2488 {
2489 b->l->v = tree->rates->cur_l[d->num];
2490 tree->rates->u_cur_l[b->num] = tree->rates->cur_l[d->num];
2491 b->l_var->v = tree->rates->cur_gamma_prior_var[d->num];
2492 }
2493
2494 if(b && (isnan(b->l->v) || isnan(b->l_var->v)))
2495 {
2496 PhyML_Fprintf(stderr,"\n. dt=%G rr=%G cr=%G ra=%G rd=%G nu=%G %f %f ",dt,rr,cr,ra,rd,nu,b->l_var->v,b->l->v);
2497 PhyML_Fprintf(stderr,"\n. ta=%G td=%G ra*cr=%G rd*cr=%G sd=%G",ta,td,ra*cr,rd*cr,SQRT(dt*nu)*cr);
2498 /* assert(FALSE); */
2499 }
2500 }
2501
2502 if(d->tax) return;
2503 else
2504 {
2505 int i;
2506 for(i=0;i<3;++i)
2507 if((d->v[i] != a) && (d->b[i] != tree->e_root))
2508 RATES_Update_Cur_Bl_Pre(d,d->v[i],d->b[i],tree);
2509 }
2510 }
2511
2512 //////////////////////////////////////////////////////////////
2513 //////////////////////////////////////////////////////////////
2514
RATES_Bl_To_Bl(t_tree * tree)2515 void RATES_Bl_To_Bl(t_tree *tree)
2516 {
2517 RATES_Bl_To_Bl_Pre(tree->n_root,tree->n_root->v[2],NULL,tree);
2518 RATES_Bl_To_Bl_Pre(tree->n_root,tree->n_root->v[1],NULL,tree);
2519 /* tree->rates->cur_l[tree->n_root->v[2]->num] = tree->a_edges[tree->e_root->num]->l->v * tree->n_root_pos; */
2520 /* tree->rates->cur_l[tree->n_root->v[1]->num] = tree->a_edges[tree->e_root->num]->l->v * (1. - tree->n_root_pos); */
2521 tree->rates->cur_l[tree->n_root->v[2]->num] = tree->a_edges[tree->e_root->num]->l->v * 0.5;
2522 tree->rates->cur_l[tree->n_root->v[1]->num] = tree->a_edges[tree->e_root->num]->l->v * (1. - 0.5);
2523 }
2524
2525 //////////////////////////////////////////////////////////////
2526 //////////////////////////////////////////////////////////////
2527
RATES_Bl_To_Bl_Pre(t_node * a,t_node * d,t_edge * b,t_tree * tree)2528 void RATES_Bl_To_Bl_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
2529 {
2530
2531 if(b)
2532 {
2533 tree->rates->cur_l[d->num] = b->l->v;
2534 }
2535
2536 if(d->tax) return;
2537 else
2538 {
2539 int i;
2540
2541 for(i=0;i<3;i++)
2542 if((d->v[i] != a) && (d->b[i] != tree->e_root))
2543 RATES_Bl_To_Bl_Pre(d,d->v[i],d->b[i],tree);
2544 }
2545 }
2546
2547 //////////////////////////////////////////////////////////////
2548 //////////////////////////////////////////////////////////////
2549
RATES_Bl_To_Ml(t_tree * tree)2550 void RATES_Bl_To_Ml(t_tree *tree)
2551 {
2552 RATES_Bl_To_Ml_Pre(tree->n_root,tree->n_root->v[2],NULL,tree);
2553 RATES_Bl_To_Ml_Pre(tree->n_root,tree->n_root->v[1],NULL,tree);
2554 tree->rates->u_ml_l[tree->e_root->num] = tree->a_edges[tree->e_root->num]->l->v;
2555 tree->rates->ml_l[tree->n_root->v[2]->num] = tree->rates->u_ml_l[tree->e_root->num] * tree->n_root_pos;
2556 tree->rates->ml_l[tree->n_root->v[1]->num] = tree->rates->u_ml_l[tree->e_root->num] * (1. - tree->n_root_pos);
2557 }
2558
2559 //////////////////////////////////////////////////////////////
2560 //////////////////////////////////////////////////////////////
2561
RATES_Bl_To_Ml_Pre(t_node * a,t_node * d,t_edge * b,t_tree * tree)2562 void RATES_Bl_To_Ml_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
2563 {
2564
2565 if(b)
2566 {
2567 tree->rates->u_ml_l[b->num] = b->l->v;
2568 tree->rates->ml_l[d->num] = b->l->v;
2569 }
2570
2571 if(d->tax) return;
2572 else
2573 {
2574 int i;
2575
2576 for(i=0;i<3;i++)
2577 if((d->v[i] != a) && (d->b[i] != tree->e_root))
2578 RATES_Bl_To_Ml_Pre(d,d->v[i],d->b[i],tree);
2579 }
2580 }
2581
2582 //////////////////////////////////////////////////////////////
2583 //////////////////////////////////////////////////////////////
2584
2585
RATES_Get_Cov_Matrix_Rooted(phydbl * unroot_cov,t_tree * tree)2586 void RATES_Get_Cov_Matrix_Rooted(phydbl *unroot_cov, t_tree *tree)
2587 {
2588 int i,dim;
2589
2590 dim = 2*tree->n_otu-3;
2591
2592 RATES_Get_Cov_Matrix_Rooted_Pre(tree->n_root,tree->n_root->v[2],NULL,unroot_cov,tree);
2593 RATES_Get_Cov_Matrix_Rooted_Pre(tree->n_root,tree->n_root->v[1],NULL,unroot_cov,tree);
2594
2595 for(i=0;i<dim+1;++i) tree->rates->cov_l[i*(dim+1)+tree->n_root->v[2]->num] /= 2.;
2596 for(i=0;i<dim+1;++i) tree->rates->cov_l[i*(dim+1)+tree->n_root->v[1]->num] /= 2.;
2597 for(i=0;i<dim+1;++i) tree->rates->cov_l[tree->n_root->v[2]->num*(dim+1)+i] /= 2.;
2598 for(i=0;i<dim+1;++i) tree->rates->cov_l[tree->n_root->v[1]->num*(dim+1)+i] /= 2.;
2599
2600 }
2601
2602 //////////////////////////////////////////////////////////////
2603 //////////////////////////////////////////////////////////////
2604
2605
RATES_Get_Cov_Matrix_Rooted_Pre(t_node * a,t_node * d,t_edge * b,phydbl * cov,t_tree * tree)2606 void RATES_Get_Cov_Matrix_Rooted_Pre(t_node *a, t_node *d, t_edge *b, phydbl *cov, t_tree *tree)
2607 {
2608 int i, dim;
2609 t_node *n;
2610
2611 dim = 2*tree->n_otu-3;
2612 n = NULL;
2613
2614 for(i=0;i<dim;i++)
2615 {
2616 if(tree->a_edges[i] != tree->e_root)
2617 {
2618 n =
2619 (tree->a_edges[i]->left->anc == tree->a_edges[i]->rght)?
2620 (tree->a_edges[i]->left):
2621 (tree->a_edges[i]->rght);
2622
2623 if(b)
2624 {
2625 tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[b->num*dim + i];
2626 }
2627 else
2628 {
2629 tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[tree->e_root->num*dim + i];
2630 }
2631 }
2632 else
2633 {
2634 n = tree->e_root->left;
2635 if(b)
2636 tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[b->num*dim + i];
2637 else
2638 tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[tree->e_root->num*dim + i];
2639
2640 n = tree->e_root->rght;
2641 if(b)
2642 tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[b->num*dim + i];
2643 else
2644 tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[tree->e_root->num*dim + i];
2645 }
2646 }
2647
2648
2649 if(d->tax) return;
2650 else
2651 {
2652 for(i=0;i<3;i++)
2653 if((d->v[i] != a) && (d->b[i] != tree->e_root))
2654 RATES_Get_Cov_Matrix_Rooted_Pre(d,d->v[i],d->b[i],cov,tree);
2655 }
2656 }
2657
2658 //////////////////////////////////////////////////////////////
2659 //////////////////////////////////////////////////////////////
2660
RATES_Covariance_Mu(t_tree * tree)2661 void RATES_Covariance_Mu(t_tree *tree)
2662 {
2663 int i,j;
2664 phydbl dt,var;
2665 int dim;
2666 int lca_num;
2667
2668 dim = 2*tree->n_otu-2;
2669
2670 For(i,dim*dim) tree->rates->cov_r[i] = 0.0;
2671
2672 dt = tree->times->nd_t[tree->n_root->v[2]->num] - tree->times->nd_t[tree->n_root->num];
2673 var = dt * tree->rates->nu;
2674 tree->rates->cov_r[tree->n_root->v[2]->num*dim+tree->n_root->v[2]->num] = var;
2675
2676
2677 dt = tree->times->nd_t[tree->n_root->v[1]->num] - tree->times->nd_t[tree->n_root->num];
2678 var = dt * tree->rates->nu;
2679 tree->rates->cov_r[tree->n_root->v[1]->num*dim+tree->n_root->v[1]->num] = var;
2680
2681 RATES_Variance_Mu_Pre(tree->n_root,tree->n_root->v[2],tree);
2682 RATES_Variance_Mu_Pre(tree->n_root,tree->n_root->v[1],tree);
2683
2684 for(i=0;i<dim;i++)
2685 {
2686 for(j=i+1;j<dim;j++)
2687 {
2688 lca_num = tree->rates->lca[i*(dim+1)+j]->num;
2689 if(lca_num < dim)
2690 {
2691 tree->rates->cov_r[i*dim+j] = tree->rates->cov_r[lca_num*dim+lca_num];
2692 tree->rates->cov_r[j*dim+i] = tree->rates->cov_r[i*dim+j];
2693 }
2694 else if(lca_num == dim)
2695 {
2696 tree->rates->cov_r[i*dim+j] = 0.0;
2697 tree->rates->cov_r[j*dim+i] = 0.0;
2698 }
2699 else
2700 {
2701 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
2702 Exit("\n");
2703 }
2704 }
2705 }
2706 }
2707
2708 //////////////////////////////////////////////////////////////
2709 //////////////////////////////////////////////////////////////
2710
2711
RATES_Variance_Mu_Pre(t_node * a,t_node * d,t_tree * tree)2712 void RATES_Variance_Mu_Pre(t_node *a, t_node *d, t_tree *tree)
2713 {
2714 int dim;
2715 phydbl var0;
2716 phydbl dt1,var1;
2717 phydbl dt2,var2;
2718 int i;
2719 int dir1, dir2;
2720
2721 dim = 2*tree->n_otu-2;
2722
2723 if(d->tax) return;
2724
2725 dir1 = dir2 = -1;
2726 for(i=0;i<3;i++)
2727 {
2728 if((d->v[i] != a) && (d->b[i] != tree->e_root))
2729 {
2730 if(dir1 < 0) dir1 = i;
2731 else dir2 = i;
2732 }
2733 }
2734
2735
2736 var0 = tree->rates->cov_r[d->num*dim+d->num];
2737
2738 dt1 = tree->times->nd_t[d->v[dir1]->num] - tree->times->nd_t[d->num];
2739 var1 = tree->rates->nu*dt1;
2740
2741 dt2 = tree->times->nd_t[d->v[dir2]->num] - tree->times->nd_t[d->num];
2742 var2 = tree->rates->nu*dt2;
2743
2744 tree->rates->cov_r[d->v[dir1]->num*dim+d->v[dir1]->num] = var0+var1;
2745 tree->rates->cov_r[d->v[dir2]->num*dim+d->v[dir2]->num] = var0+var2;
2746
2747
2748 for(i=0;i<3;i++)
2749 {
2750 if((d->v[i] != a) && (d->b[i] != tree->e_root))
2751 {
2752 RATES_Variance_Mu_Pre(d,d->v[i],tree);
2753 }
2754 }
2755 }
2756
2757 //////////////////////////////////////////////////////////////
2758 //////////////////////////////////////////////////////////////
2759
RATES_Fill_Lca_Table(t_tree * tree)2760 void RATES_Fill_Lca_Table(t_tree *tree)
2761 {
2762 int i,j;
2763 int dim;
2764
2765 dim = 2*tree->n_otu-1;
2766
2767 for(i=0;i<dim;i++)
2768 {
2769 for(j=i;j<dim;j++)
2770 {
2771 tree->rates->lca[i*dim+j] = Find_Lca_Pair_Of_Nodes(tree->a_nodes[i],tree->a_nodes[j],tree);
2772 tree->rates->lca[j*dim+i] = tree->rates->lca[i*dim+j];
2773 }
2774 }
2775 }
2776
2777 //////////////////////////////////////////////////////////////
2778 //////////////////////////////////////////////////////////////
2779
2780 /* Get V(L_{i} | L_{-i}) for all i */
RATES_Get_Conditional_Variances(t_tree * tree)2781 void RATES_Get_Conditional_Variances(t_tree *tree)
2782 {
2783 int i,j;
2784 short int *is_1;
2785 phydbl *a;
2786 int dim;
2787 t_edge *b;
2788 phydbl *cond_mu, *cond_cov;
2789
2790 dim = 2*tree->n_otu-3;
2791 a = tree->rates->_2n_vect1;
2792 is_1 = tree->rates->_2n_vect5;
2793 b = NULL;
2794 cond_mu = tree->rates->_2n_vect2;
2795 cond_cov = tree->rates->_2n2n_vect1;
2796
2797 for(i=0;i<dim;i++) a[i] = tree->rates->mean_l[i] * (Uni() * 0.2 + 0.9);
2798
2799 for(i=0;i<dim;i++)
2800 {
2801 b = tree->a_edges[i];
2802
2803 for(j=0;j<dim;j++) is_1[j] = 0;
2804 is_1[b->num] = 1;
2805
2806 For(j,dim*dim) cond_cov[j] = 0.0;
2807 for(j=0;j<dim;j++) cond_mu[j] = 0.0;
2808 Normal_Conditional(tree->rates->mean_l,tree->rates->cov_l,a,dim,is_1,1,cond_mu,cond_cov);
2809
2810 tree->rates->cond_var[b->num] = cond_cov[b->num*dim+b->num];
2811 }
2812 }
2813
2814 //////////////////////////////////////////////////////////////
2815 //////////////////////////////////////////////////////////////
2816
2817
RATES_Get_All_Reg_Coeff(t_tree * tree)2818 void RATES_Get_All_Reg_Coeff(t_tree *tree)
2819 {
2820 int i,j;
2821 short int *is_1;
2822 phydbl *a;
2823 int dim;
2824 t_edge *b;
2825
2826 dim = 2*tree->n_otu-3;
2827 a = tree->rates->_2n_vect1;
2828 is_1 = tree->rates->_2n_vect5;
2829 b = NULL;
2830
2831 for(i=0;i<dim;i++) a[i] = tree->rates->mean_l[i] * (Uni() * 0.2 + 0.9);
2832
2833 for(i=0;i<dim;i++)
2834 {
2835 b = tree->a_edges[i];
2836
2837 for(j=0;j<dim;j++) is_1[j] = 0;
2838 is_1[b->num] = 1;
2839
2840 Get_Reg_Coeff(tree->rates->mean_l,tree->rates->cov_l,a,dim,is_1,1,tree->rates->reg_coeff+b->num*dim);
2841 }
2842 }
2843
2844 //////////////////////////////////////////////////////////////
2845 //////////////////////////////////////////////////////////////
2846
2847
2848 /* Get V(L_{i} | L_{-i}) for all i */
RATES_Get_Trip_Conditional_Variances(t_tree * tree)2849 void RATES_Get_Trip_Conditional_Variances(t_tree *tree)
2850 {
2851 int i,j;
2852 short int *is_1;
2853 phydbl *a;
2854 phydbl *cond_mu, *cond_cov;
2855 t_node *n;
2856 int n_otu;
2857
2858 a = tree->rates->_2n_vect1;
2859 is_1 = tree->rates->_2n_vect5;
2860 cond_mu = tree->rates->_2n_vect2;
2861 cond_cov = tree->rates->_2n2n_vect1;
2862 n = NULL;
2863 n_otu = tree->n_otu;
2864
2865 For(i,2*n_otu-3) a[i] = tree->rates->mean_l[i] * (Uni() * 0.2 + 0.9);
2866
2867 for(i=0;i<2*n_otu-2;++i)
2868 {
2869 n = tree->a_nodes[i];
2870 if(!n->tax)
2871 {
2872 For(j,2*n_otu-3) is_1[j] = 0;
2873 is_1[n->b[0]->num] = 1;
2874 is_1[n->b[1]->num] = 1;
2875 is_1[n->b[2]->num] = 1;
2876
2877 Normal_Conditional_Unsorted(tree->rates->mean_l,tree->rates->cov_l,a,2*n_otu-3,is_1,3,cond_mu,cond_cov);
2878
2879 for(j=0;j<9;j++) tree->rates->trip_cond_cov[n->num*9+j] = cond_cov[j];
2880 }
2881 }
2882 }
2883
2884 //////////////////////////////////////////////////////////////
2885 //////////////////////////////////////////////////////////////
2886
2887
RATES_Get_All_Trip_Reg_Coeff(t_tree * tree)2888 void RATES_Get_All_Trip_Reg_Coeff(t_tree *tree)
2889 {
2890 int i,j;
2891 short int *is_1;
2892 phydbl *a;
2893 t_node *n;
2894 int n_otu;
2895
2896 a = tree->rates->_2n_vect1;
2897 is_1 = tree->rates->_2n_vect5;
2898 n = NULL;
2899 n_otu = tree->n_otu;
2900
2901 For(i,2*n_otu-3) a[i] = tree->rates->mean_l[i] * (Uni() * 0.2 + 0.9);
2902
2903 for(i=0;i<2*n_otu-2;++i)
2904 {
2905 n = tree->a_nodes[i];
2906 if(!n->tax)
2907 {
2908 For(j,2*n_otu-3) is_1[j] = 0;
2909 is_1[n->b[0]->num] = 1;
2910 is_1[n->b[1]->num] = 1;
2911 is_1[n->b[2]->num] = 1;
2912
2913 Get_Reg_Coeff(tree->rates->mean_l,tree->rates->cov_l,a,2*n_otu-3,is_1,3,tree->rates->trip_reg_coeff+n->num*(6*n_otu-9));
2914 }
2915 }
2916 }
2917
2918 //////////////////////////////////////////////////////////////
2919 //////////////////////////////////////////////////////////////
2920
RATES_Check_Lk_Rates(t_tree * tree,int * err)2921 void RATES_Check_Lk_Rates(t_tree *tree, int *err)
2922 {
2923 int i;
2924 phydbl u,u_anc;
2925 phydbl t,t_anc;
2926
2927 *err = 0;
2928
2929 For(i,2*tree->n_otu-2)
2930 {
2931 u = tree->rates->br_r[i];
2932 u_anc = tree->rates->br_r[tree->a_nodes[i]->anc->num];
2933 t = tree->times->nd_t[i];
2934 t_anc = tree->times->nd_t[tree->a_nodes[i]->anc->num];
2935
2936 if(t_anc > t)
2937 {
2938 PhyML_Printf("\n. %d %d u=%f u_anc=%f t=%f t_anc=%f",i,tree->a_nodes[i]->anc->num,u,u_anc,t,t_anc);
2939 PhyML_Printf("\n. %d %d %d",tree->n_root->num,tree->n_root->v[2]->num,tree->n_root->v[1]->num);
2940 *err = 1;
2941 }
2942 }
2943 }
2944
2945 //////////////////////////////////////////////////////////////
2946 //////////////////////////////////////////////////////////////
2947
RATES_Realized_Substitution_Rate(t_tree * tree)2948 phydbl RATES_Realized_Substitution_Rate(t_tree *tree)
2949 {
2950 return(Tree_Length(tree)/TIMES_Tree_Length(tree));
2951 }
2952
2953 //////////////////////////////////////////////////////////////
2954 //////////////////////////////////////////////////////////////
2955
RATES_Expected_Tree_Length(t_tree * tree)2956 phydbl RATES_Expected_Tree_Length(t_tree *tree)
2957 {
2958 int n;
2959 phydbl mean;
2960
2961 n = 0;
2962 mean = 0.0;
2963 RATES_Expected_Tree_Length_Pre(tree->n_root,tree->n_root->v[2],1.0,&mean,&n,tree);
2964 RATES_Expected_Tree_Length_Pre(tree->n_root,tree->n_root->v[1],1.0,&mean,&n,tree);
2965
2966 if(n != 2*tree->n_otu-2)
2967 {
2968 PhyML_Fprintf(stderr,"\n. n=%d 2n-2=%d",n,2*tree->n_otu-2);
2969 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
2970 Exit("\n");
2971 }
2972
2973 return mean;
2974 }
2975
2976 //////////////////////////////////////////////////////////////
2977 //////////////////////////////////////////////////////////////
2978
RATES_Expected_Tree_Length_Pre(t_node * a,t_node * d,phydbl eranc,phydbl * mean,int * n,t_tree * tree)2979 void RATES_Expected_Tree_Length_Pre(t_node *a, t_node *d, phydbl eranc, phydbl *mean, int *n, t_tree *tree)
2980 {
2981 phydbl erdes;
2982 int i;
2983 phydbl loc_mean;
2984 int loc_n;
2985
2986
2987 /* erdes = u_anc - */
2988 /* sd*(Dnorm((tree->rates->min_rate-u_anc)/sd,.0,1.) - Dnorm((tree->rates->max_rate-u_anc)/sd,.0,1.))/ */
2989 /* (Pnorm((tree->rates->max_rate-u_anc)/sd,.0,1.) - Pnorm((tree->rates->min_rate-u_anc)/sd,.0,1.)); */
2990
2991 /* erdes = Norm_Trunc_Mean(eranc,sd,tree->rates->min_rate,tree->rates->max_rate); */
2992 erdes = 1.0;
2993
2994 loc_mean = *mean;
2995 loc_n = *n;
2996
2997 loc_mean *= loc_n;
2998 loc_mean += erdes;
2999 loc_mean /= (loc_n + 1);
3000
3001 *mean = loc_mean;
3002 *n = loc_n + 1;
3003
3004
3005 if(d->tax) return;
3006 else
3007 for(i=0;i<3;i++)
3008 if(d->v[i] != a && d->b[i] != tree->e_root)
3009 RATES_Expected_Tree_Length_Pre(d,d->v[i],erdes,mean,n,tree);
3010 }
3011
3012 //////////////////////////////////////////////////////////////
3013 //////////////////////////////////////////////////////////////
3014
RATES_Update_Norm_Fact(t_tree * tree)3015 void RATES_Update_Norm_Fact(t_tree *tree)
3016 {
3017 /* int i; */
3018 /* phydbl r,t,t_anc; */
3019 /* phydbl num,denom; */
3020
3021 /* num = denom = 0.0; */
3022
3023 /* For(i,2*tree->n_otu-2) */
3024 /* { */
3025 /* r = tree->rates->br_r[i]; */
3026 /* t = tree->times->nd_t[i]; */
3027 /* t_anc = tree->times->nd_t[tree->a_nodes[i]->anc->num]; */
3028
3029 /* num += (t-t_anc); */
3030 /* denom += (t-t_anc) * r; */
3031 /* } */
3032 /* tree->rates->norm_fact = num/denom; */
3033
3034 tree->rates->norm_fact = 1.0;
3035 }
3036
3037 //////////////////////////////////////////////////////////////
3038 //////////////////////////////////////////////////////////////
3039
RATES_Normalise_Rates(t_tree * tree)3040 void RATES_Normalise_Rates(t_tree *tree)
3041 {
3042 phydbl expr,curr;
3043 int i;
3044
3045 curr = RATES_Average_Substitution_Rate(tree);
3046 curr /= tree->rates->clock_r;
3047
3048 /* Set expected mean rate to one such that clock_r is
3049 easy to interpret regarding the mean values of br_r */
3050 expr = 1.0;
3051
3052 For(i,2*tree->n_otu-2)
3053 {
3054 tree->rates->br_r[i] *= expr/curr;
3055
3056 if(tree->rates->br_r[i] > tree->rates->max_rate)
3057 tree->rates->br_r[i] = tree->rates->max_rate;
3058
3059 if(tree->rates->br_r[i] < tree->rates->min_rate)
3060 tree->rates->br_r[i] = tree->rates->min_rate;
3061 }
3062
3063 tree->rates->clock_r *= curr/expr;
3064 /* Branch lengths therefore do not change */
3065
3066 /* if(tree->rates->clock_r < tree->rates->min_clock) */
3067 /* { */
3068 /* PhyML_Printf("\n. Curr mean rates: %G",curr); */
3069 /* PhyML_Printf("\n. Set clock rate to: %G",tree->rates->clock_r); */
3070 /* tree->rates->clock_r = tree->rates->min_clock; */
3071 /* } */
3072 /* if(tree->rates->clock_r > tree->rates->max_clock) */
3073 /* { */
3074 /* PhyML_Printf("\n. Curr mean rates: %G",curr); */
3075 /* PhyML_Printf("\n. Set clock rate to: %G",tree->rates->clock_r); */
3076 /* tree->rates->clock_r = tree->rates->max_clock; */
3077 /* } */
3078
3079 RATES_Update_Norm_Fact(tree);
3080 RATES_Update_Cur_Bl(tree);
3081 }
3082
3083 //////////////////////////////////////////////////////////////
3084 //////////////////////////////////////////////////////////////
3085
3086
RATES_Find_Max_Dt_Bisec(phydbl r,phydbl r_mean,phydbl ta,phydbl tc,phydbl nu,phydbl threshp,int inf)3087 phydbl RATES_Find_Max_Dt_Bisec(phydbl r, phydbl r_mean, phydbl ta, phydbl tc, phydbl nu, phydbl threshp, int inf)
3088 {
3089 phydbl cdfr,cdf0;
3090 phydbl sd;
3091 phydbl trunc_cdf;
3092 phydbl ori_tc, ori_ta;
3093 phydbl tb;
3094
3095 ori_tc = tc;
3096 ori_ta = ta;
3097
3098 /* PhyML_Printf("\n Max %s r=%f r_mean%f",inf?"inf":"sup",r,r_mean); */
3099 do
3100 {
3101 tb = ta + (tc - ta)/2.;
3102
3103
3104 sd = SQRT(nu*(ori_tc - tb));
3105 cdfr = Pnorm(r ,r_mean,sd);
3106 cdf0 = Pnorm(.0,r_mean,sd);
3107
3108 if(inf)
3109 trunc_cdf = (cdfr - cdf0)/(1. - cdf0);
3110 else
3111 trunc_cdf = (1. - cdfr)/(1. - cdf0);
3112
3113 /* PhyML_Printf("\n. ta=%15f tb=%15f tc=%15f cdf = %15f",ta,tb,tc,trunc_cdf); */
3114
3115 if(trunc_cdf > threshp)
3116 {
3117 ta = tb;
3118 }
3119 else
3120 {
3121 tc = tb;
3122 }
3123
3124 }while((tc - ta)/(ori_tc - ori_ta) > 0.001);
3125
3126 if(tb < ori_ta)
3127 {
3128 PhyML_Fprintf(stderr,"\n. tb < ta r=%f r_mean=%f",r,r_mean);
3129 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
3130 Exit("\n");
3131 }
3132 if(tb > ori_tc)
3133 {
3134 PhyML_Fprintf(stderr,"\n. tb > tc r=%f r_mean=%f ori_ta=%f ori_tc=%f tb=%f",r,r_mean,ori_ta,ori_tc,tb);
3135 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
3136 Exit("\n");
3137 }
3138
3139 return tb;
3140 }
3141
3142 //////////////////////////////////////////////////////////////
3143 //////////////////////////////////////////////////////////////
3144
3145
RATES_Find_Min_Dt_Bisec(phydbl r,phydbl r_mean,phydbl ta,phydbl tc,phydbl nu,phydbl threshp,int inf)3146 phydbl RATES_Find_Min_Dt_Bisec(phydbl r, phydbl r_mean, phydbl ta, phydbl tc, phydbl nu, phydbl threshp, int inf)
3147 {
3148 phydbl cdfr,cdf0;
3149 phydbl sd;
3150 phydbl trunc_cdf;
3151 phydbl ori_tc, ori_ta;
3152 phydbl tb;
3153
3154 ori_tc = tc;
3155 ori_ta = ta;
3156
3157 /* PhyML_Printf("\n Min %s r=%f r_mean=%f",inf?"inf":"sup",r,r_mean); */
3158 do
3159 {
3160
3161 tb = ta + (tc - ta)/2.;
3162
3163
3164 sd = SQRT(nu*(tb - ori_ta));
3165 cdfr = Pnorm(r ,r_mean,sd);
3166 cdf0 = Pnorm(.0,r_mean,sd);
3167
3168 if(inf)
3169 trunc_cdf = (cdfr - cdf0)/(1. - cdf0);
3170 else
3171 trunc_cdf = (1. - cdfr)/(1. - cdf0);
3172
3173 /* PhyML_Printf("\n. ta=%15f tb=%15f tc=%15f cdf = %15f",ta,tb,tc,trunc_cdf); */
3174
3175 if(trunc_cdf > threshp)
3176 {
3177 tc = tb;
3178 }
3179 else
3180 {
3181 ta = tb;
3182 }
3183
3184 }while((tc - ta)/(ori_tc - ori_ta) > 0.001);
3185
3186 if(tb < ori_ta)
3187 {
3188 PhyML_Fprintf(stderr,"\n. tb < ta r=%f r_mean=%f",r,r_mean);
3189 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
3190 Exit("\n");
3191 }
3192 if(tb > ori_tc)
3193 {
3194 PhyML_Fprintf(stderr,"\n. tb > tc r=%f r_mean=%f ori_ta=%f ori_tc=%f tb=%f",r,r_mean,ori_ta,ori_tc,tb);
3195 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
3196 Exit("\n");
3197 }
3198
3199 return tb;
3200 }
3201
3202 //////////////////////////////////////////////////////////////
3203 //////////////////////////////////////////////////////////////
3204
3205
RATES_Min_Max_Interval(phydbl u0,phydbl u1,phydbl u2,phydbl u3,phydbl t0,phydbl t2,phydbl t3,phydbl * t_min,phydbl * t_max,phydbl nu,phydbl p_thresh,t_tree * tree)3206 void RATES_Min_Max_Interval(phydbl u0, phydbl u1, phydbl u2, phydbl u3, phydbl t0, phydbl t2, phydbl t3,
3207 phydbl *t_min, phydbl *t_max, phydbl nu, phydbl p_thresh, t_tree *tree)
3208 {
3209 int n_eval;
3210 int i;
3211 phydbl sum_cdf;
3212 phydbl *cdf;
3213 phydbl t1;
3214 phydbl mint2t3;
3215
3216 mint2t3 = MIN(t2,t3);
3217 n_eval = 5;
3218
3219 cdf = (phydbl *)mCalloc(n_eval,sizeof(phydbl));
3220
3221 sum_cdf = .0;
3222 for(i=0;i<n_eval;i++)
3223 {
3224 t1 = t0 + (i + 1)*(mint2t3 - t0)/(n_eval + 1);
3225 cdf[i] =
3226 Dnorm_Trunc(u1,u0,SQRT(nu*(t1-t0)),tree->rates->min_rate,tree->rates->max_rate) *
3227 Dnorm_Trunc(u2,u1,SQRT(nu*(t2-t1)),tree->rates->min_rate,tree->rates->max_rate) *
3228 Dnorm_Trunc(u3,u1,SQRT(nu*(t3-t1)),tree->rates->min_rate,tree->rates->max_rate) *
3229 (mint2t3 - t0) / (n_eval + 1);
3230 if(i) cdf[i] += cdf[i-1];
3231 }
3232 sum_cdf = cdf[i-1];
3233
3234 for(i=0;i<n_eval;i++) cdf[i] /= sum_cdf;
3235
3236 /* PhyML_Printf("\n"); */
3237 /* for(i=0;i<n_eval;i++) PhyML_Printf("\n* %f %f",cdf[i],sum_cdf); */
3238
3239 for(i=0;i<n_eval;i++)
3240 if(cdf[i] > p_thresh)
3241 {
3242 *t_min = t0 + (i + 1)*(mint2t3 - t0)/(n_eval + 1);
3243 break;
3244 }
3245
3246 for(i=0;i<n_eval;i++)
3247 if(cdf[i] > (1. - p_thresh))
3248 {
3249 *t_max = t0 + (i + 1)*(mint2t3 - t0)/(n_eval + 1);
3250 break;
3251 }
3252
3253 Free(cdf);
3254 }
3255
3256 //////////////////////////////////////////////////////////////
3257 //////////////////////////////////////////////////////////////
3258
3259
RATES_Get_Correction_Factor(phydbl mode,phydbl sd,int * err,t_tree * tree)3260 phydbl RATES_Get_Correction_Factor(phydbl mode, phydbl sd, int *err, t_tree *tree)
3261 {
3262
3263 phydbl K,X0,X1,X2,Y0,Y1,Y2;
3264 phydbl eps=0.01;
3265 phydbl A,B;
3266 phydbl slope,inter;
3267 phydbl denom;
3268
3269 *err = NO;
3270
3271
3272 /* DO NOTHING */
3273 return 0.0;
3274
3275
3276
3277 A = tree->rates->min_rate / sd;
3278 B = tree->rates->max_rate / sd;
3279 K = mode / sd;
3280
3281 X0 = 0.;
3282 /* Y0 = Dnorm(X0-K,.0,1.)/(1. - Pnorm(X0-K,.0,1.)) - X0; */
3283 Y0 = (Dnorm(A-K+X0,.0,1.)-Dnorm(B-K+X0,.0,1.))/(Pnorm(B-K+X0,.0,1.)-Pnorm(A-K+X0,.0,1.)) - X0;
3284
3285 X1 = .1;
3286 /* Y1 = Dnorm(X1-K,.0,1.)/(1. - Pnorm(X1-K,.0,1.)) - X1; */
3287 Y1 = (Dnorm(A-K+X1,.0,1.)-Dnorm(B-K+X1,.0,1.))/(Pnorm(B-K+X1,.0,1.)-Pnorm(A-K+X1,.0,1.)) - X1;
3288
3289 /* printf("\n. ^^ mean=%f sd=%f",mode,sd); */
3290
3291 /* printf("\n. X0=%f Y0=%f X1=%f Y1=%f",X0,Y0,X1,Y1); */
3292
3293 do
3294 {
3295
3296 slope = (Y1-Y0)/(X1-X0);
3297 inter = Y0 - X0*(Y1-Y0)/(X1-X0);
3298
3299 X2 = -inter/slope;
3300
3301 denom = (Pnorm(B-K+X2,.0,1.)-Pnorm(A-K+X2,.0,1.));
3302
3303 if(denom < 1.E-10)
3304 {
3305 /* printf("\n. X2 = %f Y2=%f num=%f denom=%f Y1=%f Y0=%f X1=%f X0=%f mode=%f sd=%f",X2,Y2,num,denom,Y1,Y0,X1,X0,mode,sd); */
3306 *err = YES;
3307 break;
3308 }
3309
3310 /* Y2 = Dnorm(X2-K,.0,1.)/(1. - Pnorm(X2-K,.0,1.)) - X2; */
3311 Y2 = (Dnorm(A-K+X2,.0,1.)-Dnorm(B-K+X2,.0,1.))/(Pnorm(B-K+X2,.0,1.)-Pnorm(A-K+X2,.0,1.)) - X2;
3312
3313 /* printf("\n. X2 = %f Y2=%f num=%f denom=%f Y1=%f Y0=%f X1=%f X0=%f",X2,Y2,num,denom,Y1,Y0,X1,X0); */
3314
3315 if(X2 > X1)
3316 {
3317 X0 = X1;
3318 X1 = X2;
3319 Y0 = Y1;
3320 Y1 = Y2;
3321 }
3322 else
3323 {
3324 X1 = X0;
3325 X0 = X2;
3326 Y1 = Y0;
3327 Y0 = Y2;
3328 }
3329 }while(fabs(Y2) > eps);
3330
3331 /* printf("\n. shift = %f X2=%f Y2 = %f",X2*sd,X2,Y2); */
3332
3333 return X2 * sd;
3334 }
3335
3336 //////////////////////////////////////////////////////////////
3337 //////////////////////////////////////////////////////////////
3338
3339
Sample_Average_Rate(t_node * a,t_node * d,t_tree * tree)3340 phydbl Sample_Average_Rate(t_node *a, t_node *d, t_tree *tree)
3341 {
3342 return(-1.);
3343 }
3344
3345 //////////////////////////////////////////////////////////////
3346 //////////////////////////////////////////////////////////////
3347
RATES_Update_Mean_Br_Len(int iter,t_tree * tree)3348 void RATES_Update_Mean_Br_Len(int iter, t_tree *tree)
3349 {
3350 int i,dim;
3351 phydbl *mean;
3352
3353 if(tree->rates->update_mean_l == NO) return;
3354
3355 dim = 2*tree->n_otu-3;
3356 mean = tree->rates->mean_l;
3357
3358 for(i=0;i<dim;i++)
3359 {
3360 mean[i] *= (phydbl)iter;
3361 mean[i] += tree->a_edges[i]->l->v;
3362 mean[i] /= (phydbl)(iter+1);
3363 }
3364 }
3365
3366 //////////////////////////////////////////////////////////////
3367 //////////////////////////////////////////////////////////////
3368
3369
RATES_Update_Cov_Br_Len(int iter,t_tree * tree)3370 void RATES_Update_Cov_Br_Len(int iter, t_tree *tree)
3371 {
3372 int i,j,dim;
3373 phydbl *mean,*cov;
3374
3375 if(tree->rates->update_cov_l == NO) return;
3376
3377 dim = 2*tree->n_otu-3;
3378 mean = tree->rates->mean_l;
3379 cov = tree->rates->cov_l;
3380
3381 for(i=0;i<dim;i++)
3382 {
3383 for(j=0;j<dim;j++)
3384 {
3385 cov[i*dim+j] += mean[i]*mean[j];
3386 cov[i*dim+j] *= (phydbl)tree->mcmc->run;
3387 cov[i*dim+j] += tree->a_edges[i]->l->v * tree->a_edges[j]->l->v;
3388 cov[i*dim+j] /= (phydbl)(tree->mcmc->run+1);
3389 cov[i*dim+j] -= mean[i]*mean[j];
3390
3391 if(i == j && cov[i*dim+j] < MIN_VAR_BL) cov[i*dim+j] = MIN_VAR_BL;
3392 }
3393 }
3394 }
3395
3396 //////////////////////////////////////////////////////////////
3397 //////////////////////////////////////////////////////////////
3398
RATES_Set_Mean_L(t_tree * tree)3399 void RATES_Set_Mean_L(t_tree *tree)
3400 {
3401 int i;
3402 For(i,2*tree->n_otu-3)
3403 {
3404 tree->rates->mean_l[i] = tree->a_edges[i]->l->v;
3405 }
3406 }
3407
3408 //////////////////////////////////////////////////////////////
3409 //////////////////////////////////////////////////////////////
3410
RATES_Record_Times(t_tree * mixt_tree)3411 void RATES_Record_Times(t_tree *mixt_tree)
3412 {
3413 int i;
3414 t_tree *tree;
3415
3416 tree = mixt_tree;
3417 do
3418 {
3419 if(tree->times->nd_t_recorded == YES)
3420 {
3421 PhyML_Fprintf(stderr,"\n. Overwriting recorded times is forbidden.\n");
3422 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
3423 Exit("\n");
3424 }
3425
3426 For(i,2*tree->n_otu-1) tree->times->buff_t[i] = tree->times->nd_t[i];
3427 tree = tree->next;
3428 }
3429 while(tree);
3430 }
3431
3432 //////////////////////////////////////////////////////////////
3433 //////////////////////////////////////////////////////////////
3434
3435
RATES_Reset_Times(t_tree * mixt_tree)3436 void RATES_Reset_Times(t_tree *mixt_tree)
3437 {
3438 int i;
3439 t_tree *tree;
3440
3441 tree = mixt_tree;
3442 do
3443 {
3444 tree->times->nd_t_recorded = NO;
3445 for(i=0;i<2*tree->n_otu-1;++i) tree->times->nd_t[i] = tree->times->buff_t[i];
3446 tree = tree->next;
3447 }
3448 while(tree);
3449 }
3450
3451 //////////////////////////////////////////////////////////////
3452 //////////////////////////////////////////////////////////////
3453
RATES_Record_Rates(t_tree * tree)3454 void RATES_Record_Rates(t_tree *tree)
3455 {
3456 int i;
3457
3458 if(tree->rates->br_r_recorded == YES)
3459 {
3460 PhyML_Fprintf(stderr,"\n. Overwriting recorded rates is forbidden.\n");
3461 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
3462 Exit("\n");
3463 }
3464
3465 for(i=0;i<2*tree->n_otu-2;++i) tree->rates->buff_br_r[i] = tree->rates->br_r[i];
3466 for(i=0;i<2*tree->n_otu-1;++i) tree->rates->buff_nd_r[i] = tree->rates->nd_r[i];
3467 }
3468
3469 //////////////////////////////////////////////////////////////
3470 //////////////////////////////////////////////////////////////
3471
3472
RATES_Reset_Rates(t_tree * tree)3473 void RATES_Reset_Rates(t_tree *tree)
3474 {
3475 int i;
3476 tree->rates->br_r_recorded = NO;
3477 for(i=0;i<2*tree->n_otu-2;++i) tree->rates->br_r[i] = tree->rates->buff_br_r[i];
3478 for(i=0;i<2*tree->n_otu-1;++i) tree->rates->nd_r[i] = tree->rates->buff_nd_r[i];
3479 }
3480
3481 //////////////////////////////////////////////////////////////
3482 //////////////////////////////////////////////////////////////
3483
RATES_Set_Clock_And_Nu_Max(t_tree * tree)3484 void RATES_Set_Clock_And_Nu_Max(t_tree *tree)
3485 {
3486 phydbl dt,nu;
3487 phydbl min_t;
3488 int i;
3489 phydbl step;
3490 phydbl l_max;
3491 phydbl max_clock;
3492 phydbl r_max;
3493 phydbl tune;
3494 phydbl pa,pb;
3495
3496 if(tree->rates->model == THORNE ||
3497 tree->rates->model == GUINDON)
3498 {
3499 tune = 1.05;
3500 r_max = tree->rates->max_rate;
3501 l_max = tree->mod->l_max;
3502
3503 min_t = .0;
3504 for(i=0;i<2*tree->n_otu-1;++i) if(tree->times->t_prior_min[i] < min_t) min_t = tree->times->t_prior_min[i];
3505
3506 dt = FABS(min_t);
3507 max_clock = l_max / dt;
3508
3509 nu = 1.E-10;
3510 step = 1.E-1;
3511 do
3512 {
3513 do
3514 {
3515 nu += step;
3516 pa = Dnorm(0.0, 0.0,SQRT(nu*dt));
3517 pb = Dnorm(r_max,0.0,SQRT(nu*dt));
3518 }while(pa/pb > tune);
3519 nu -= step;
3520 step /= 10.;
3521 }while(step > 1.E-10);
3522
3523 tree->rates->max_nu = nu;
3524 /* tree->rates->max_nu = 1.0; */
3525 tree->rates->max_clock = max_clock;
3526
3527 PhyML_Printf("\n. Clock rate parameter upper bound set to %f expected subst./site/time unit",tree->rates->max_clock);
3528 PhyML_Printf("\n. Autocorrelation parameter upper bound set to %f",tree->rates->max_nu);
3529 }
3530 }
3531
3532 //////////////////////////////////////////////////////////////
3533 //////////////////////////////////////////////////////////////
3534
RATES_Set_Birth_Rate_Boundaries(t_tree * tree)3535 void RATES_Set_Birth_Rate_Boundaries(t_tree *tree)
3536 {
3537 phydbl lbda;
3538 phydbl p_above_min,p_below_max;
3539 phydbl min,max;
3540 int assign = YES;
3541
3542 min = -0.01*tree->times->t_prior_max[tree->n_root->num];
3543 max = -100.*tree->times->t_prior_min[tree->n_root->num];
3544
3545 for(lbda = 0.0001; lbda < 10; lbda+=0.0001)
3546 {
3547 p_above_min = 1. - POW(1.-exp(-lbda*min),tree->n_otu);
3548 p_below_max = POW(1.-exp(-lbda*max),tree->n_otu);
3549
3550 if(p_above_min < 1.E-10)
3551 {
3552 tree->times->birth_rate_max = lbda;
3553 break;
3554 }
3555 if(p_below_max > 1.E-10 && assign==YES)
3556 {
3557 assign = NO;
3558 tree->times->birth_rate_min = lbda;
3559 }
3560 }
3561
3562 /* tree->times->birth_rate_min = 1.E-6; */
3563 /* tree->times->birth_rate_max = 1.; */
3564 PhyML_Printf("\n. Birth rate lower bound set to %f.",tree->times->birth_rate_min);
3565 PhyML_Printf("\n. Birth rate upper bound set to %f.",tree->times->birth_rate_max);
3566
3567 }
3568
3569
3570 //////////////////////////////////////////////////////////////
3571 //////////////////////////////////////////////////////////////
3572
RATES_Get_Mean_Rate_In_Subtree(t_node * root,t_tree * tree)3573 phydbl RATES_Get_Mean_Rate_In_Subtree(t_node *root, t_tree *tree)
3574 {
3575 phydbl sum;
3576 int n;
3577
3578 sum = 0.0;
3579 n = 0;
3580
3581 if(root->tax == NO)
3582 {
3583 if(root == tree->n_root)
3584 {
3585 RATES_Get_Mean_Rate_In_Subtree_Pre(root,root->v[2],&sum,&n,tree);
3586 RATES_Get_Mean_Rate_In_Subtree_Pre(root,root->v[1],&sum,&n,tree);
3587 }
3588 else
3589 {
3590 int i;
3591 for(i=0;i<3;i++)
3592 {
3593 if(root->v[i] != root->anc && root->b[i] != tree->e_root)
3594 {
3595 RATES_Get_Mean_Rate_In_Subtree_Pre(root,root->v[i],&sum,&n,tree);
3596 }
3597 }
3598 }
3599 return sum/(phydbl)n;
3600 }
3601 else
3602 {
3603 return 0.0;
3604 }
3605 }
3606
3607 //////////////////////////////////////////////////////////////
3608 //////////////////////////////////////////////////////////////
3609
3610
RATES_Get_Mean_Rate_In_Subtree_Pre(t_node * a,t_node * d,phydbl * sum,int * n,t_tree * tree)3611 void RATES_Get_Mean_Rate_In_Subtree_Pre(t_node *a, t_node *d, phydbl *sum, int *n, t_tree *tree)
3612 {
3613 /* (*sum) += exp(tree->rates->nd_r[d->num]); */
3614
3615 if(tree->rates->model == LOGNORMAL ||
3616 tree->rates->model == THORNE ||
3617 tree->rates->model == STRICTCLOCK)
3618 {
3619 (*sum) += tree->rates->br_r[d->num];
3620 }
3621 else if(tree->rates->model == GUINDON)
3622 {
3623 (*sum) += tree->rates->nd_r[d->num];
3624 }
3625
3626 else assert(FALSE);
3627
3628 (*n) += 1;
3629
3630 if(d->tax == YES) return;
3631 else
3632 {
3633 int i;
3634 for(i=0;i<3;++i)
3635 {
3636 if(d->v[i] != a && d->b[i] != tree->e_root)
3637 {
3638 RATES_Get_Mean_Rate_In_Subtree_Pre(d,d->v[i],sum,n,tree);
3639 }
3640 }
3641 }
3642 }
3643
3644 //////////////////////////////////////////////////////////////
3645 //////////////////////////////////////////////////////////////
3646
RATES_Get_Model_Name(int model)3647 char *RATES_Get_Model_Name(int model)
3648 {
3649 char *s;
3650
3651 s = (char *)mCalloc(T_MAX_NAME,sizeof(char));
3652
3653 switch(model)
3654 {
3655 case GUINDON : {strcpy(s,"geometric Brownian"); break;}
3656 case THORNE : {strcpy(s,"gbd"); break;}
3657 case LOGNORMAL : {strcpy(s,"lognormal"); break;}
3658 case STRICTCLOCK : {strcpy(s,"strict clock"); break;}
3659 default :
3660 {
3661 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
3662 Exit("\n");
3663 break;
3664 }
3665 }
3666 return s;
3667 }
3668
3669 //////////////////////////////////////////////////////////////
3670 //////////////////////////////////////////////////////////////
3671
RATES_Check_Edge_Length_Consistency(t_tree * mixt_tree)3672 int RATES_Check_Edge_Length_Consistency(t_tree *mixt_tree)
3673 {
3674 int i;
3675 t_tree *tree;
3676
3677 tree = mixt_tree;
3678
3679 do
3680 {
3681 if(tree->is_mixt_tree == YES) tree = tree->next;
3682
3683 for(i=0;i<2*tree->n_otu-3;++i)
3684 {
3685 if(tree->a_edges[i]->left->anc == tree->a_edges[i]->rght)
3686 {
3687 if(Are_Equal(tree->rates->cur_l[tree->a_edges[i]->left->num],
3688 tree->a_edges[i]->l->v,
3689 1.E-5) == NO)
3690 {
3691 PhyML_Fprintf(stderr,"\n. cur_l: %G l: %G is_root: %d",
3692 tree->rates->cur_l[tree->a_edges[i]->left->num],
3693 tree->a_edges[i]->l->v,
3694 tree->a_edges[i] == tree->e_root);
3695 return 0;
3696 }
3697 }
3698
3699 if(tree->a_edges[i]->rght->anc == tree->a_edges[i]->left)
3700 {
3701 if(Are_Equal(tree->rates->cur_l[tree->a_edges[i]->rght->num],
3702 tree->a_edges[i]->l->v,
3703 1.E-5) == NO)
3704 {
3705 PhyML_Fprintf(stderr,"\n. cur_l: %G l: %G is_root: %d",
3706 tree->rates->cur_l[tree->a_edges[i]->rght->num],
3707 tree->a_edges[i]->l->v,
3708 tree->a_edges[i] == tree->e_root);
3709 return 0;
3710 }
3711 }
3712 }
3713
3714 tree = tree->next;
3715 }
3716 while(tree);
3717
3718 return 1;
3719
3720 }
3721
3722 //////////////////////////////////////////////////////////////
3723 //////////////////////////////////////////////////////////////
3724
3725 //////////////////////////////////////////////////////////////
3726 //////////////////////////////////////////////////////////////
3727
3728 //////////////////////////////////////////////////////////////
3729 //////////////////////////////////////////////////////////////
3730
3731 //////////////////////////////////////////////////////////////
3732 //////////////////////////////////////////////////////////////
3733
3734
3735
3736