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
14 /* Routines for molecular clock trees and molecular dating */
15
16
17 #include "times.h"
18
19 //////////////////////////////////////////////////////////////
20 //////////////////////////////////////////////////////////////
21
TIMES_Least_Square_Node_Times(t_edge * e_root,t_tree * tree)22 void TIMES_Least_Square_Node_Times(t_edge *e_root, t_tree *tree)
23 {
24
25 /* Solve A.x = b, where x is the vector of times estimated
26 under the least square criterion and b is the vector
27 of edge lengths
28
29 A is a n x n matrix, with n being the number of
30 nodes in a rooted tree (i.e. 2*n_otu-1).
31
32 */
33
34 phydbl *A, *b, *x;
35 int n;
36 int i,j;
37 t_node *root;
38
39 n = 2*tree->n_otu-1;
40
41 A = (phydbl *)mCalloc(n*n,sizeof(phydbl));
42 b = (phydbl *)mCalloc(n, sizeof(phydbl));
43 x = (phydbl *)mCalloc(n, sizeof(phydbl));
44
45 if(!tree->n_root && e_root) Add_Root(e_root,tree);
46 else if(!e_root) Add_Root(tree->a_edges[0],tree);
47
48 root = tree->n_root;
49
50 TIMES_Least_Square_Node_Times_Pre(root,root->v[1],A,b,n,tree);
51 TIMES_Least_Square_Node_Times_Pre(root,root->v[2],A,b,n,tree);
52
53 b[root->num] = tree->e_root->l->v/2.;
54
55 A[root->num * n + root->num] = 1.0;
56 A[root->num * n + root->v[2]->num] = -.5;
57 A[root->num * n + root->v[1]->num] = -.5;
58
59 if(!Matinv(A, n, n,YES))
60 {
61 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s').\n",__FILE__,__LINE__,__FUNCTION__);
62 Exit("\n");
63 }
64
65 for(i=0;i<n;i++) x[i] = .0;
66 for(i=0;i<n;i++) for(j=0;j<n;j++) x[i] += A[i*n+j] * b[j];
67
68 for(i=0;i<n-1;i++) tree->times->nd_t[tree->a_nodes[i]->num] = -x[i];
69 tree->times->nd_t[root->num] = -x[n-1];
70 tree->n_root->b[2]->l->v = tree->times->nd_t[root->v[2]->num] - tree->times->nd_t[root->num];
71 tree->n_root->b[1]->l->v = tree->times->nd_t[root->v[1]->num] - tree->times->nd_t[root->num];
72
73 Free(A);
74 Free(b);
75 Free(x);
76
77 return;
78
79 }
80
81 //////////////////////////////////////////////////////////////
82 //////////////////////////////////////////////////////////////
83
TIMES_Least_Square_Node_Times_Pre(t_node * a,t_node * d,phydbl * A,phydbl * b,int n,t_tree * tree)84 void TIMES_Least_Square_Node_Times_Pre(t_node *a, t_node *d, phydbl *A, phydbl *b, int n, t_tree *tree)
85 {
86 if(d->tax)
87 {
88 A[d->num * n + d->num] = 1.;
89
90 /* Set the time stamp at tip nodes to 0.0 */
91 /* PhyML_Printf("\n. Tip t_node date set to 0"); */
92 b[d->num] = 0.0;
93 return;
94 }
95 else
96 {
97 int i;
98
99 for(i=0;i<3;++i)
100 if((d->v[i] != a) && (d->b[i] != tree->e_root))
101 TIMES_Least_Square_Node_Times_Pre(d,d->v[i],A,b,n,tree);
102
103 A[d->num * n + d->num] = 1.;
104 b[d->num] = .0;
105 for(i=0;i<3;++i)
106 {
107 A[d->num * n + d->v[i]->num] = -1./3.;
108 if(d->v[i] != a) b[d->num] += d->b[i]->l->v;
109 else b[d->num] -= d->b[i]->l->v;
110 }
111 b[d->num] /= 3.;
112 }
113 }
114
115 //////////////////////////////////////////////////////////////
116 //////////////////////////////////////////////////////////////
117
118
119 /* Adjust t_node times in order to have correct time stamp ranking with
120 respect to the tree topology */
121
TIMES_Adjust_Node_Times(t_tree * tree)122 void TIMES_Adjust_Node_Times(t_tree *tree)
123 {
124 TIMES_Adjust_Node_Times_Pre(tree->n_root->v[2],tree->n_root->v[1],tree);
125 TIMES_Adjust_Node_Times_Pre(tree->n_root->v[1],tree->n_root->v[2],tree);
126
127 if(tree->times->nd_t[tree->n_root->num] > MIN(tree->times->nd_t[tree->n_root->v[2]->num],
128 tree->times->nd_t[tree->n_root->v[1]->num]))
129 {
130 tree->times->nd_t[tree->n_root->num] = MIN(tree->times->nd_t[tree->n_root->v[2]->num],
131 tree->times->nd_t[tree->n_root->v[1]->num]);
132 }
133 }
134
135 //////////////////////////////////////////////////////////////
136 //////////////////////////////////////////////////////////////
137
138
TIMES_Adjust_Node_Times_Pre(t_node * a,t_node * d,t_tree * tree)139 void TIMES_Adjust_Node_Times_Pre(t_node *a, t_node *d, t_tree *tree)
140 {
141 if(d->tax) return;
142 else
143 {
144 int i;
145 phydbl min_height;
146
147 for(i=0;i<3;i++)
148 if((d->v[i] != a) && (d->b[i] != tree->e_root))
149 {
150 TIMES_Adjust_Node_Times_Pre(d,d->v[i],tree);
151 }
152
153 min_height = 0.0;
154 for(i=0;i<3;i++)
155 {
156 if((d->v[i] != a) && (d->b[i] != tree->e_root))
157 {
158 if(tree->times->nd_t[d->v[i]->num] < min_height)
159 {
160 min_height = tree->times->nd_t[d->v[i]->num];
161 }
162 }
163 }
164
165 if(tree->times->nd_t[d->num] > min_height) tree->times->nd_t[d->num] = min_height;
166
167 if(tree->times->nd_t[d->num] < -100.) tree->times->nd_t[d->num] = -100.;
168
169 }
170 }
171
172 //////////////////////////////////////////////////////////////
173 //////////////////////////////////////////////////////////////
174
175
176 /* Multiply each time stamp at each internal
177 t_node by 'tree->time_stamp_mult'.
178 */
179
TIMES_Mult_Time_Stamps(t_tree * tree)180 void TIMES_Mult_Time_Stamps(t_tree *tree)
181 {
182 int i;
183 For(i,2*tree->n_otu-2) tree->times->nd_t[tree->a_nodes[i]->num] *= FABS(tree->mod->s_opt->tree_size_mult);
184 tree->times->nd_t[tree->n_root->num] *= FABS(tree->mod->s_opt->tree_size_mult);
185 }
186
187 //////////////////////////////////////////////////////////////
188 //////////////////////////////////////////////////////////////
189
TIMES_Print_Node_Times(t_node * a,t_node * d,t_tree * tree)190 void TIMES_Print_Node_Times(t_node *a, t_node *d, t_tree *tree)
191 {
192 t_edge *b;
193 int i;
194
195 b = NULL;
196 for(i=0;i<3;i++) if((d->v[i]) && (d->v[i] == a)) {b = d->b[i]; break;}
197
198 PhyML_Printf("\n. (%3d %3d) a->t = %12f d->t = %12f (#=%12f) b->l->v = %12f [%12f;%12f]",
199 a->num,d->num,
200 tree->times->nd_t[a->num],
201 tree->times->nd_t[d->num],
202 tree->times->nd_t[a->num]-tree->times->nd_t[d->num],
203 (b)?(b->l->v):(-1.0),
204 tree->times->t_prior_min[d->num],
205 tree->times->t_prior_max[d->num]);
206 if(d->tax) return;
207 else
208 {
209 int i;
210 for(i=0;i<3;i++)
211 if((d->v[i] != a) && (d->b[i] != tree->e_root))
212 TIMES_Print_Node_Times(d,d->v[i],tree);
213 }
214 }
215
216 //////////////////////////////////////////////////////////////
217 //////////////////////////////////////////////////////////////
218
TIMES_Set_All_Node_Priors(t_tree * tree)219 void TIMES_Set_All_Node_Priors(t_tree *tree)
220 {
221 int i;
222 phydbl min_prior;
223
224 /* Set all t_prior_max values */
225 TIMES_Set_All_Node_Priors_Bottom_Up(tree->n_root,tree->n_root->v[2],tree);
226 TIMES_Set_All_Node_Priors_Bottom_Up(tree->n_root,tree->n_root->v[1],tree);
227
228 tree->times->t_prior_max[tree->n_root->num] =
229 MIN(tree->times->t_prior_max[tree->n_root->num],
230 MIN(tree->times->t_prior_max[tree->n_root->v[2]->num],
231 tree->times->t_prior_max[tree->n_root->v[1]->num]));
232
233
234 /* Set all t_prior_min values */
235 if(!tree->times->t_has_prior[tree->n_root->num])
236 {
237 min_prior = 1.E+10;
238 for(i=0;i<2*tree->n_otu-2;++i)
239 {
240 if(tree->times->t_has_prior[i])
241 {
242 if(tree->times->t_prior_min[i] < min_prior)
243 min_prior = tree->times->t_prior_min[i];
244 }
245 }
246 tree->times->t_prior_min[tree->n_root->num] = 2.0 * min_prior;
247 /* tree->times->t_prior_min[tree->n_root->num] = 10. * min_prior; */
248 }
249
250 if(tree->times->t_prior_min[tree->n_root->num] > 0.0)
251 {
252 PhyML_Fprintf(stderr,"\n. Failed to set the lower bound for the root node.");
253 PhyML_Fprintf(stderr,"\n. Make sure at least one of the calibration interval");
254 PhyML_Fprintf(stderr,"\n. provides a lower bound.");
255 Exit("\n");
256 }
257
258
259 TIMES_Set_All_Node_Priors_Top_Down(tree->n_root,tree->n_root->v[2],tree);
260 TIMES_Set_All_Node_Priors_Top_Down(tree->n_root,tree->n_root->v[1],tree);
261
262 Get_Node_Ranks(tree);
263 TIMES_Set_Floor(tree);
264 }
265
266 //////////////////////////////////////////////////////////////
267 //////////////////////////////////////////////////////////////
268
269
TIMES_Set_All_Node_Priors_Bottom_Up(t_node * a,t_node * d,t_tree * tree)270 void TIMES_Set_All_Node_Priors_Bottom_Up(t_node *a, t_node *d, t_tree *tree)
271 {
272 int i;
273 phydbl t_sup;
274
275 if(d->tax) return;
276 else
277 {
278 t_node *v1, *v2; /* the two sons of d */
279
280 for(i=0;i<3;i++)
281 {
282 if((d->v[i] != a) && (d->b[i] != tree->e_root))
283 {
284 TIMES_Set_All_Node_Priors_Bottom_Up(d,d->v[i],tree);
285 }
286 }
287
288 v1 = v2 = NULL;
289 for(i=0;i<3;i++) if((d->v[i] != a) && (d->b[i] != tree->e_root))
290 {
291 if(!v1) v1 = d->v[i];
292 else v2 = d->v[i];
293 }
294
295 if(tree->times->t_has_prior[d->num] == YES)
296 {
297 t_sup = MIN(tree->times->t_prior_max[d->num],
298 MIN(tree->times->t_prior_max[v1->num],
299 tree->times->t_prior_max[v2->num]));
300
301 tree->times->t_prior_max[d->num] = t_sup;
302
303 if(tree->times->t_prior_max[d->num] < tree->times->t_prior_min[d->num])
304 {
305 PhyML_Fprintf(stderr,"\n. prior_min=%f prior_max=%f",tree->times->t_prior_min[d->num],tree->times->t_prior_max[d->num]);
306 PhyML_Fprintf(stderr,"\n. Inconsistency in the prior settings detected at node %d",d->num);
307 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function %s)\n\n",__FILE__,__LINE__,__FUNCTION__);
308 Warn_And_Exit("\n");
309 }
310 }
311 else
312 {
313 tree->times->t_prior_max[d->num] =
314 MIN(tree->times->t_prior_max[v1->num],
315 tree->times->t_prior_max[v2->num]);
316 }
317 }
318 }
319
320 //////////////////////////////////////////////////////////////
321 //////////////////////////////////////////////////////////////
322
323
TIMES_Set_All_Node_Priors_Top_Down(t_node * a,t_node * d,t_tree * tree)324 void TIMES_Set_All_Node_Priors_Top_Down(t_node *a, t_node *d, t_tree *tree)
325 {
326 if(d->tax) return;
327 else
328 {
329 int i;
330
331 if(tree->times->t_has_prior[d->num] == YES)
332 {
333 tree->times->t_prior_min[d->num] = MAX(tree->times->t_prior_min[d->num],tree->times->t_prior_min[a->num]);
334
335 if(tree->times->t_prior_max[d->num] < tree->times->t_prior_min[d->num])
336 {
337 PhyML_Fprintf(stderr,"\n. prior_min=%f prior_max=%f",tree->times->t_prior_min[d->num],tree->times->t_prior_max[d->num]);
338 PhyML_Fprintf(stderr,"\n. Inconsistency in the prior settings detected at t_node %d",d->num);
339 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function %s)\n\n",__FILE__,__LINE__,__FUNCTION__);
340 Warn_And_Exit("\n");
341 }
342 }
343 else
344 {
345 tree->times->t_prior_min[d->num] = tree->times->t_prior_min[a->num];
346 }
347
348 for(i=0;i<3;i++)
349 {
350 if((d->v[i] != a) && (d->b[i] != tree->e_root))
351 {
352 TIMES_Set_All_Node_Priors_Top_Down(d,d->v[i],tree);
353 }
354 }
355 }
356 }
357
358 //////////////////////////////////////////////////////////////
359 //////////////////////////////////////////////////////////////
360
TIMES_Set_Floor(t_tree * tree)361 void TIMES_Set_Floor(t_tree *tree)
362 {
363 TIMES_Set_Floor_Post(tree->n_root,tree->n_root->v[2],tree);
364 TIMES_Set_Floor_Post(tree->n_root,tree->n_root->v[1],tree);
365 tree->times->t_floor[tree->n_root->num] = MIN(tree->times->t_floor[tree->n_root->v[2]->num],
366 tree->times->t_floor[tree->n_root->v[1]->num]);
367 }
368
369 //////////////////////////////////////////////////////////////
370 //////////////////////////////////////////////////////////////
371
TIMES_Set_Floor_Post(t_node * a,t_node * d,t_tree * tree)372 void TIMES_Set_Floor_Post(t_node *a, t_node *d, t_tree *tree)
373 {
374 if(d->tax)
375 {
376 tree->times->t_floor[d->num] = tree->times->nd_t[d->num];
377 d->rank_max = d->rank;
378 return;
379 }
380 else
381 {
382 int i;
383 t_node *v1,*v2;
384
385 v1 = v2 = NULL;
386 for(i=0;i<3;i++)
387 {
388 if(d->v[i] != a && d->b[i] != tree->e_root)
389 {
390 TIMES_Set_Floor_Post(d,d->v[i],tree);
391 if(!v1) v1 = d->v[i];
392 else v2 = d->v[i];
393 }
394 }
395 tree->times->t_floor[d->num] = MIN(tree->times->t_floor[v1->num],
396 tree->times->t_floor[v2->num]);
397
398 if(tree->times->t_floor[v1->num] < tree->times->t_floor[v2->num])
399 {
400 d->rank_max = v1->rank_max;
401 }
402 else if(tree->times->t_floor[v2->num] < tree->times->t_floor[v1->num])
403 {
404 d->rank_max = v2->rank_max;
405 }
406 else
407 {
408 d->rank_max = MAX(v1->rank_max,v2->rank_max);
409 }
410 }
411 }
412
413 //////////////////////////////////////////////////////////////
414 //////////////////////////////////////////////////////////////
415
416 /* Does it work for serial samples? */
TIMES_Log_Conditional_Uniform_Density(t_tree * tree)417 phydbl TIMES_Log_Conditional_Uniform_Density(t_tree *tree)
418 {
419 phydbl min,max;
420 phydbl dens;
421 int i;
422
423 min = tree->times->nd_t[tree->n_root->num];
424
425 dens = 0.0;
426 For(i,2*tree->n_otu-1)
427 {
428 if((tree->a_nodes[i]->tax == NO) && (tree->a_nodes[i] != tree->n_root))
429 {
430 max = tree->times->t_floor[i];
431
432 dens += log(Dorder_Unif(tree->times->nd_t[i],
433 tree->a_nodes[i]->rank-1,
434 tree->a_nodes[i]->rank_max-2,
435 min,max));
436 }
437 }
438 return dens;
439 }
440
441 //////////////////////////////////////////////////////////////
442 //////////////////////////////////////////////////////////////
443 // Returns the marginal density of tree height assuming the
444 // Yule model of speciation.
TIMES_Lk_Yule_Root_Marginal(t_tree * tree)445 phydbl TIMES_Lk_Yule_Root_Marginal(t_tree *tree)
446 {
447 int n;
448 int j;
449 t_node *nd;
450 phydbl *t,*ts;
451 phydbl lbda;
452 phydbl T;
453
454 lbda = tree->times->birth_rate;
455 t = tree->times->nd_t;
456 ts = tree->times->time_slice_lims;
457 T = ts[0] - t[tree->n_root->num];
458
459 n = 0;
460 nd = NULL;
461 For(j,2*tree->n_otu-2)
462 {
463 nd = tree->a_nodes[j];
464
465 if((t[nd->num] > ts[0] && t[nd->anc->num] < ts[0]) || // lineage that is crossing ts[0]
466 (nd->tax == YES && Are_Equal(t[nd->num],ts[0],1.E-6) == YES)) // tip that is lying on ts[0]
467 n++;
468 }
469
470 return LnGamma(n+1) + log(lbda) - 2.*lbda*T + (n-2.)*log(1. - exp(-lbda*T));
471 }
472
473 //////////////////////////////////////////////////////////////
474 //////////////////////////////////////////////////////////////
475 // Returns the joint density of internal node heights assuming
476 // the Yule model of speciation.
TIMES_Lk_Yule_Joint(t_tree * tree)477 phydbl TIMES_Lk_Yule_Joint(t_tree *tree)
478 {
479 int i,j;
480 phydbl loglk;
481 phydbl *t;
482 phydbl dt;
483 int n; // number of lineages at a given time point
484 phydbl lbda;
485 t_node *nd;
486 phydbl *ts;
487 int *tr;
488 phydbl top_t;
489 short int *interrupted;
490 phydbl sumdt;
491
492 interrupted = (short int *)mCalloc(tree->n_otu,sizeof(short int));
493
494 t = tree->times->nd_t;
495 ts = tree->times->time_slice_lims;
496 tr = tree->times->t_rank;
497 lbda = tree->times->birth_rate;
498
499 TIMES_Update_Node_Ordering(tree);
500
501 for(j=0;j<tree->n_otu;j++) interrupted[j] = NO;
502
503 loglk = .0;
504
505 sumdt = .0;
506 n = 1;
507 For(i,2*tree->n_otu-2) // t[tr[0]] is the time of the oldest node, t[tr[1]], the second oldest and so on...
508 {
509
510 for(j=0;j<tree->n_otu;j++)
511 if((t[j] < t[tr[i]]) && (interrupted[j] == NO))
512 {
513 interrupted[j] = YES;
514 n--; // How many lineages have stopped above t[tr[i]]?
515 }
516
517 top_t = t[tr[i+1]];
518 dt = top_t - t[tr[i]];
519 sumdt += dt;
520
521 /* printf("\n. %d node up=%d [%f] node do=%d [%f] dt=%f",i,tr[i],t[tr[i]],tr[i+1],t[tr[i+1]],dt); */
522
523 if(n<1)
524 {
525 PhyML_Fprintf(stderr,"\n. i=%d tr[i]=%f",i,t[tr[i]]);
526 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
527 Exit("\n");
528 }
529
530 if(dt > 1.E-10) loglk += log((n+1)*lbda) - (n+1)*lbda*dt;
531 n++;
532 }
533
534 /* printf("\n. sumdt = %f th=%f",sumdt,tree->times->nd_t[tree->n_root->num]); */
535 /* printf("\n0 loglk = %f",loglk); */
536
537 for(i=0;i<tree->times->n_time_slices-1;i++)
538 {
539 n = 0;
540 dt = 0.;
541 For(j,2*tree->n_otu-2)
542 {
543 nd = tree->a_nodes[j];
544 if(t[nd->num] > ts[i] && t[nd->anc->num] < ts[i]) // How many lineages are crossing this time slice limit?
545 {
546 n++;
547 if(t[nd->num] < dt) dt = t[nd->num]; // take the oldest node younger than the time slice
548 }
549 }
550 dt -= ts[i];
551 loglk += log(n*lbda) - n*lbda*dt;
552 }
553
554 /* printf("\n1 loglk = %f",loglk); */
555
556 Free(interrupted);
557
558 return loglk;
559 }
560
561 //////////////////////////////////////////////////////////////
562 //////////////////////////////////////////////////////////////
563 // Returns the conditional density of internal node heights
564 // given the tree height under the Yule model. Uses the order
565 // statistics 'simplification' as described in Yang and Rannala, 2005.
TIMES_Lk_Yule_Order(t_tree * tree)566 phydbl TIMES_Lk_Yule_Order(t_tree *tree)
567 {
568 int j;
569 phydbl *t,*tf;
570 t_node *n;
571 phydbl loglk;
572 phydbl loglbda;
573 phydbl lbda;
574 phydbl *tp_min,*tp_max;
575 phydbl lower_bound,upper_bound;
576 /* phydbl root_height; */
577
578 tp_min = tree->times->t_prior_min;
579 tp_max = tree->times->t_prior_max;
580 tf = tree->times->t_floor;
581 t = tree->times->nd_t;
582 n = NULL;
583 loglbda = log(tree->times->birth_rate);
584 lbda = tree->times->birth_rate;
585 lower_bound = -1.;
586 upper_bound = -1.;
587 /* root_height = FABS(tree->times->nd_t[tree->n_root->num]); */
588
589 /*! Adapted from Equation (6) in T. Stadler's Systematic Biology, 2012 paper with
590 sampling fraction set to 1 and death rate set to 0. Dropped the 1/(n-1) scaling
591 factor. */
592
593 /* loglk = 0.0; */
594 /* For(j,2*tree->n_otu-2) */
595 /* { */
596 /* n = tree->a_nodes[j]; */
597 /* lower_bound = MAX(FABS(tf[j]),FABS(tp_max[j])); */
598 /* upper_bound = MIN(FABS(t[tree->n_root->num]),FABS(tp_min[j])); */
599
600 /* if(n->tax == NO) */
601 /* { */
602 /* loglk += (loglbda - lbda * FABS(t[j])); */
603 /* /\* loglk -= log(exp(-lbda*lower_bound) - exp(-lbda*upper_bound)); // incorporate calibration boundaries here. *\/ */
604 /* } */
605 /* } */
606
607
608 /*! Adapted from Equation (7) in T. Stadler's Systematic Biology, 2012 paper with
609 sampling fraction set to 1 and death rate set to 0. */
610
611 // Check that each node is within its calibration-derived interval
612 For(j,2*tree->n_otu-1) if(t[j] < tp_min[j] || t[j] > tp_max[j]) return(-INFINITY);
613
614 loglk = 0.0;
615 For(j,2*tree->n_otu-2)
616 {
617 n = tree->a_nodes[j];
618 lower_bound = MAX(FABS(tf[j]),FABS(tp_max[j]));
619 upper_bound = FABS(tp_min[j]);
620
621 if(n->tax == NO)
622 {
623 loglk += (loglbda - lbda * FABS(t[j]));
624 loglk -= log(exp(-lbda*lower_bound) - exp(-lbda*upper_bound)); // incorporate calibration boundaries here.
625 }
626
627 if(isinf(loglk) || isnan(loglk))
628 {
629 /* PhyML_Printf("\n. Lower bound: %f",lower_bound); */
630 /* PhyML_Printf("\n. Upper bound: %f",upper_bound); */
631 /* PhyML_Printf("\n. tf: %f tp_max: %f tp_min: %f ",tf[j],tp_max[j],tp_min[j]); */
632 /* PhyML_Printf("\n. exp1: %f",exp(-lbda*lower_bound)); */
633 /* PhyML_Printf("\n. exp2: %f",exp(-lbda*upper_bound)); */
634 /* PhyML_Printf("\n. diff: %f",exp(-lbda*lower_bound) - exp(-lbda*upper_bound)); */
635 /* Exit("\n"); */
636 return(-INFINITY);
637 }
638 }
639
640 lower_bound = MAX(FABS(tf[tree->n_root->num]),FABS(tp_max[tree->n_root->num]));
641 upper_bound = FABS(tp_min[tree->n_root->num]);
642 loglk += log(2) + loglbda - 2.*lbda * FABS(t[tree->n_root->num]);
643 loglk -= log(exp(-2.*lbda*lower_bound) - exp(-2.*lbda*upper_bound));
644
645 if(isinf(loglk) || isnan(loglk))
646 {
647 /* PhyML_Printf("\n. * Lower bound: %f",lower_bound); */
648 /* PhyML_Printf("\n. * Upper bound: %f",upper_bound); */
649 /* PhyML_Printf("\n. * tf: %f tp_max: %f tp_min: %f",tf[tree->n_root->num],tp_max[tree->n_root->num],tp_min[tree->n_root->num]); */
650 /* PhyML_Printf("\n. * exp1: %f",exp(-2.*lbda*lower_bound)); */
651 /* PhyML_Printf("\n. * exp2: %f",exp(-2.*lbda*upper_bound)); */
652 /* PhyML_Printf("\n. * diff: %f",exp(-2.*lbda*lower_bound) - exp(-2.*lbda*upper_bound)); */
653 /* Exit("\n"); */
654 return(-INFINITY);
655 }
656
657
658 return(loglk);
659 }
660
661 //////////////////////////////////////////////////////////////
662 //////////////////////////////////////////////////////////////
663
TIMES_Lk_Times(int verbose,t_tree * tree)664 phydbl TIMES_Lk_Times(int verbose, t_tree *tree)
665 {
666 DATE_Assign_Primary_Calibration(tree);
667 DATE_Update_T_Prior_MinMax(tree);
668 tree->times->c_lnL_times = TIMES_Lk_Birth_Death(verbose,tree);
669 return(tree->times->c_lnL_times);
670 }
671
672 //////////////////////////////////////////////////////////////
673 //////////////////////////////////////////////////////////////
674
675
TIMES_Lk_Times_Trav(t_node * a,t_node * d,phydbl lim_inf,phydbl lim_sup,phydbl * logdens,t_tree * tree)676 void TIMES_Lk_Times_Trav(t_node *a, t_node *d, phydbl lim_inf, phydbl lim_sup, phydbl *logdens, t_tree *tree)
677 {
678 int i;
679
680 if(!d->tax)
681 {
682 /* if(tree->times->nd_t[d->num] > lim_sup) */
683 /* { */
684 /* lim_inf = lim_sup; */
685 /* lim_sup = 0.0; */
686 /* For(i,2*tree->n_otu-2) */
687 /* if((tree->times->t_floor[i] < lim_sup) && (tree->times->t_floor[i] > tree->times->nd_t[d->num])) */
688 /* lim_sup = tree->times->t_floor[i]; */
689 /* } */
690
691 /* if(tree->times->nd_t[d->num] < lim_inf || tree->times->nd_t[d->num] > lim_sup) */
692 /* { */
693 /* PhyML_Printf("\n. nd_t = %f lim_inf = %f lim_sup = %f",tree->times->nd_t[d->num],lim_inf,lim_sup); */
694 /* PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__); */
695 /* Exit("\n"); */
696 /* } */
697
698 lim_inf = tree->times->nd_t[tree->n_root->num];
699 lim_sup = tree->times->t_floor[d->num];
700
701 *logdens = *logdens + log(lim_sup - lim_inf);
702 }
703
704 if(d->tax == YES) return;
705 else
706 {
707 for(i=0;i<3;i++)
708 {
709 if(d->v[i] != a && d->b[i] != tree->e_root)
710 {
711 TIMES_Lk_Times_Trav(d,d->v[i],lim_inf,lim_sup,logdens,tree);
712 }
713 }
714 }
715 }
716
717 //////////////////////////////////////////////////////////////
718 //////////////////////////////////////////////////////////////
719
TIMES_Log_Number_Of_Ranked_Labelled_Histories(t_node * root,int per_slice,t_tree * tree)720 phydbl TIMES_Log_Number_Of_Ranked_Labelled_Histories(t_node *root, int per_slice, t_tree *tree)
721 {
722 int i;
723 phydbl logn;
724 t_node *v1,*v2;
725 int n1,n2;
726
727 TIMES_Update_Curr_Slice(tree);
728
729 logn = .0;
730 v1 = v2 = NULL;
731 if(root == tree->n_root)
732 {
733 TIMES_Log_Number_Of_Ranked_Labelled_Histories_Post(root,root->v[2],per_slice,&logn,tree);
734 TIMES_Log_Number_Of_Ranked_Labelled_Histories_Post(root,root->v[1],per_slice,&logn,tree);
735 v1 = root->v[2];
736 v2 = root->v[1];
737 }
738 else
739 {
740 for(i=0;i<3;i++)
741 {
742 if(root->v[i] != root->anc && root->b[i] != tree->e_root)
743 {
744 TIMES_Log_Number_Of_Ranked_Labelled_Histories_Post(root,root->v[i],per_slice,&logn,tree);
745 if(!v1) v1 = root->v[i];
746 else v2 = root->v[i];
747 }
748 }
749 }
750
751
752 if(per_slice == NO)
753 {
754 n1 = tree->rates->n_tips_below[v1->num];
755 n2 = tree->rates->n_tips_below[v2->num];
756 }
757 else
758 {
759 if(tree->times->curr_slice[v1->num] == tree->times->curr_slice[root->num])
760 n1 = tree->rates->n_tips_below[v1->num];
761 else
762 n1 = 1;
763
764 if(tree->times->curr_slice[v2->num] == tree->times->curr_slice[root->num])
765 n2 = tree->rates->n_tips_below[v2->num];
766 else
767 n2 = 1;
768 }
769
770 tree->rates->n_tips_below[root->num] = n1+n2;
771
772 logn += Factln(n1+n2-2) - Factln(n1-1) - Factln(n2-1);
773
774 return(logn);
775 }
776
777 //////////////////////////////////////////////////////////////
778 //////////////////////////////////////////////////////////////
779
780
TIMES_Log_Number_Of_Ranked_Labelled_Histories_Post(t_node * a,t_node * d,int per_slice,phydbl * logn,t_tree * tree)781 void TIMES_Log_Number_Of_Ranked_Labelled_Histories_Post(t_node *a, t_node *d, int per_slice, phydbl *logn, t_tree *tree)
782 {
783 if(d->tax == YES)
784 {
785 tree->rates->n_tips_below[d->num] = 1;
786 return;
787 }
788 else
789 {
790 int i,n1,n2;
791 t_node *v1, *v2;
792
793 for(i=0;i<3;i++)
794 {
795 if(d->v[i] != a && d->b[i] != tree->e_root)
796 {
797 TIMES_Log_Number_Of_Ranked_Labelled_Histories_Post(d,d->v[i],per_slice,logn,tree);
798 }
799 }
800
801 v1 = v2 = NULL;
802 for(i=0;i<3;i++)
803 {
804 if(d->v[i] != a && d->b[i] != tree->e_root)
805 {
806 if(v1 == NULL) {v1 = d->v[i];}
807 else {v2 = d->v[i];}
808 }
809 }
810
811
812 if(per_slice == NO)
813 {
814 n1 = tree->rates->n_tips_below[v1->num];
815 n2 = tree->rates->n_tips_below[v2->num];
816 }
817 else
818 {
819 if(tree->times->curr_slice[v1->num] == tree->times->curr_slice[d->num])
820 n1 = tree->rates->n_tips_below[v1->num];
821 else
822 n1 = 1;
823
824 if(tree->times->curr_slice[v2->num] == tree->times->curr_slice[d->num])
825 n2 = tree->rates->n_tips_below[v2->num];
826 else
827 n2 = 1;
828 }
829
830 tree->rates->n_tips_below[d->num] = n1+n2;
831
832 (*logn) += Factln(n1+n2-2) - Factln(n1-1) - Factln(n2-1);
833 }
834 }
835
836 //////////////////////////////////////////////////////////////
837 //////////////////////////////////////////////////////////////
838
839
TIMES_Update_Curr_Slice(t_tree * tree)840 void TIMES_Update_Curr_Slice(t_tree *tree)
841 {
842 int i,j;
843
844 for(i=0;i<2*tree->n_otu-1;++i)
845 {
846 for(j=0;j<tree->times->n_time_slices;j++)
847 {
848 if(!(tree->times->nd_t[i] > tree->times->time_slice_lims[j])) break;
849 }
850 tree->times->curr_slice[i] = j;
851
852 /* PhyML_Printf("\n. Node %3d [%12f] is in slice %3d.",i,tree->times->nd_t[i],j); */
853 }
854 }
855
856 //////////////////////////////////////////////////////////////
857 //////////////////////////////////////////////////////////////
858
TIMES_Wrap_Lk_Coalescent(t_edge * b,t_tree * tree,supert_tree * stree)859 phydbl TIMES_Wrap_Lk_Coalescent(t_edge *b, t_tree *tree, supert_tree *stree)
860 {
861 return TIMES_Lk_Coalescent(tree);
862 }
863
864 /*////////////////////////////////////////////////////////////
865 ////////////////////////////////////////////////////////////*/
866
867
TIMES_Lk_Coalescent(t_tree * tree)868 phydbl TIMES_Lk_Coalescent(t_tree *tree)
869 {
870 t_node *n;
871 int n_lineages;
872 phydbl lnP;
873
874 Get_Node_Ranks_From_Times(tree);
875
876
877 lnP = 0.0;
878 n_lineages = 1;
879 n = tree->n_root;
880 while(n->rk_next)
881 {
882 if(n->tax == YES) n_lineages--;
883 else n_lineages++;
884
885 lnP += -n_lineages * (n_lineages-1.) / (2.*tree->times->scaled_pop_size) * fabs(tree->times->nd_t[n->num] - tree->times->nd_t[n->rk_next->num]);
886
887 n = n->rk_next;
888 }
889
890 lnP -= (tree->n_otu-1) * log(tree->times->scaled_pop_size);
891
892 tree->times->c_lnL_times = lnP;
893
894 return(lnP);
895
896 }
897
898 //////////////////////////////////////////////////////////////
899 //////////////////////////////////////////////////////////////
900
TIMES_Lk_Uniform_Core(t_tree * tree)901 phydbl TIMES_Lk_Uniform_Core(t_tree *tree)
902 {
903 phydbl logn;
904
905 logn = TIMES_Log_Number_Of_Ranked_Labelled_Histories(tree->n_root,YES,tree);
906
907 tree->times->c_lnL_times = 0.0;
908 TIMES_Lk_Uniform_Post(tree->n_root,tree->n_root->v[2],tree);
909 TIMES_Lk_Uniform_Post(tree->n_root,tree->n_root->v[1],tree);
910
911 /* printf("\n. ^ %f %f %f", */
912 /* (phydbl)(tree->rates->n_tips_below[tree->n_root->num]-2.), */
913 /* log(tree->times->time_slice_lims[tree->times->curr_slice[tree->n_root->num]] - */
914 /* tree->times->nd_t[tree->n_root->num]), */
915 /* (phydbl)(tree->rates->n_tips_below[tree->n_root->num]-2.) * */
916 /* log(tree->times->time_slice_lims[tree->times->curr_slice[tree->n_root->num]] - */
917 /* tree->times->nd_t[tree->n_root->num])); */
918
919 tree->times->c_lnL_times +=
920 Factln(tree->rates->n_tips_below[tree->n_root->num]-2.) -
921 (phydbl)(tree->rates->n_tips_below[tree->n_root->num]-2.) *
922 log(tree->times->time_slice_lims[tree->times->curr_slice[tree->n_root->num]] -
923 tree->times->nd_t[tree->n_root->num]);
924
925 tree->times->c_lnL_times -= logn;
926
927 return(tree->times->c_lnL_times);
928 }
929
930 //////////////////////////////////////////////////////////////
931 //////////////////////////////////////////////////////////////
932
933
TIMES_Get_Number_Of_Time_Slices(t_tree * tree)934 void TIMES_Get_Number_Of_Time_Slices(t_tree *tree)
935 {
936 int i;
937
938 tree->times->n_time_slices=0;
939 TIMES_Get_Number_Of_Time_Slices_Post(tree->n_root,tree->n_root->v[2],tree);
940 TIMES_Get_Number_Of_Time_Slices_Post(tree->n_root,tree->n_root->v[1],tree);
941 Qksort(tree->times->time_slice_lims,NULL,0,tree->times->n_time_slices-1);
942
943 if(tree->times->n_time_slices > 1)
944 {
945 PhyML_Printf("\n");
946 PhyML_Printf("\n. Sequences were collected at %d different time points.",tree->times->n_time_slices);
947 for(i=0;i<tree->times->n_time_slices;i++) printf("\n+ [%3d] time point @ %12f ",i+1,tree->times->time_slice_lims[i]);
948 }
949 }
950
951 //////////////////////////////////////////////////////////////
952 //////////////////////////////////////////////////////////////
953
954
TIMES_Get_Number_Of_Time_Slices_Post(t_node * a,t_node * d,t_tree * tree)955 void TIMES_Get_Number_Of_Time_Slices_Post(t_node *a, t_node *d, t_tree *tree)
956 {
957 int i;
958
959 if(d->tax == YES)
960 {
961 for(i=0;i<tree->times->n_time_slices;i++)
962 if(Are_Equal(tree->times->t_floor[d->num],tree->times->time_slice_lims[i],1.E-6)) break;
963
964 if(i == tree->times->n_time_slices)
965 {
966 tree->times->time_slice_lims[i] = tree->times->t_floor[d->num];
967 tree->times->n_time_slices++;
968 }
969 return;
970 }
971 else
972 {
973 for(i=0;i<3;i++)
974 if(d->v[i] != a && d->b[i] != tree->e_root)
975 TIMES_Get_Number_Of_Time_Slices_Post(d,d->v[i],tree);
976 }
977 }
978
979 //////////////////////////////////////////////////////////////
980 //////////////////////////////////////////////////////////////
981
982
TIMES_Get_N_Slice_Spans(t_tree * tree)983 void TIMES_Get_N_Slice_Spans(t_tree *tree)
984 {
985 int i,j;
986
987 For(i,2*tree->n_otu-2)
988 {
989 if(tree->a_nodes[i]->tax == NO)
990 {
991 for(j=0;j<tree->times->n_time_slices;j++)
992 {
993 if(Are_Equal(tree->times->t_floor[i],tree->times->time_slice_lims[j],1.E-6))
994 {
995 tree->times->n_time_slice_spans[i] = j+1;
996 /* PhyML_Printf("\n. Node %3d spans %3d slices [%12f].", */
997 /* i+1, */
998 /* tree->rates->n_slice_spans[i], */
999 /* tree->times->t_floor[i]); */
1000 break;
1001 }
1002 }
1003 }
1004 }
1005 }
1006
1007 //////////////////////////////////////////////////////////////
1008 //////////////////////////////////////////////////////////////
1009
1010
TIMES_Lk_Uniform_Post(t_node * a,t_node * d,t_tree * tree)1011 void TIMES_Lk_Uniform_Post(t_node *a, t_node *d, t_tree *tree)
1012 {
1013 if(d->tax == YES) return;
1014 else
1015 {
1016 int i;
1017
1018 for(i=0;i<3;i++)
1019 {
1020 if(d->v[i] != a && d->b[i] != tree->e_root)
1021 {
1022 TIMES_Lk_Uniform_Post(d,d->v[i],tree);
1023 }
1024 }
1025
1026 if(tree->times->curr_slice[a->num] != tree->times->curr_slice[d->num])
1027 {
1028 tree->times->c_lnL_times +=
1029 Factln(tree->rates->n_tips_below[d->num]-1.) -
1030 (phydbl)(tree->rates->n_tips_below[d->num]-1.) *
1031 log(tree->times->time_slice_lims[tree->times->curr_slice[d->num]] -
1032 tree->times->nd_t[d->num]);
1033 }
1034 }
1035 }
1036
1037 //////////////////////////////////////////////////////////////
1038 //////////////////////////////////////////////////////////////
1039
1040 /* Set the root position so that most of the taxa in the outgroup
1041 correspond to the most ancient time point.
1042 */
TIMES_Set_Root_Given_Tip_Dates(t_tree * tree)1043 void TIMES_Set_Root_Given_Tip_Dates(t_tree *tree)
1044 {
1045 int i,j;
1046 t_node *left,*rght;
1047 int n_left_in, n_left_out;
1048 int n_rght_in, n_rght_out;
1049 t_edge *b,*best;
1050 phydbl eps,score,max_score;
1051
1052 Free_Bip(tree);
1053 Alloc_Bip(tree);
1054 Get_Bip(tree->a_nodes[0],tree->a_nodes[0]->v[0],tree);
1055
1056 left = rght = NULL;
1057 b = best = NULL;
1058 n_left_in = n_left_out = -1;
1059 n_rght_in = n_rght_out = -1;
1060 eps = 1.E-6;
1061 score = max_score = -1.;
1062
1063 For(i,2*tree->n_otu-3)
1064 {
1065 left = tree->a_edges[i]->left;
1066 rght = tree->a_edges[i]->rght;
1067 b = tree->a_edges[i];
1068
1069 n_left_in = 0;
1070 For(j,left->bip_size[b->l_r])
1071 if(FABS(tree->times->nd_t[left->bip_node[b->l_r][j]->num] - tree->times->time_slice_lims[0]) < eps)
1072 n_left_in++;
1073
1074 n_left_out = left->bip_size[b->l_r]-n_left_in;
1075
1076 n_rght_in = 0;
1077 For(j,rght->bip_size[b->r_l])
1078 if(FABS(tree->times->nd_t[rght->bip_node[b->r_l][j]->num] - tree->times->time_slice_lims[0]) < eps)
1079 n_rght_in++;
1080
1081 n_rght_out = rght->bip_size[b->r_l]-n_rght_in;
1082
1083
1084 /* score = POW((phydbl)(n_left_in)/(phydbl)(n_left_in+n_left_out)- */
1085 /* (phydbl)(n_rght_in)/(phydbl)(n_rght_in+n_rght_out),2); */
1086 /* score = (phydbl)(n_left_in * n_rght_out + eps)/(n_left_out * n_rght_in + eps); */
1087 /* score = (phydbl)(n_left_in * n_rght_out + eps); */
1088 score = FABS((phydbl)((n_left_in+1.) * (n_rght_out+1.)) - (phydbl)((n_left_out+1.) * (n_rght_in+1.)));
1089
1090 if(score > max_score)
1091 {
1092 max_score = score;
1093 best = b;
1094 }
1095 }
1096
1097 Add_Root(best,tree);
1098 }
1099
1100 //////////////////////////////////////////////////////////////
1101 //////////////////////////////////////////////////////////////
1102
1103 /* Update the ranking of node heights. Use bubble sort algorithm */
1104 /* t_rank[i] is the node number that has rank i */
1105
TIMES_Update_Node_Ordering(t_tree * tree)1106 void TIMES_Update_Node_Ordering(t_tree *tree)
1107 {
1108 int buff;
1109 int i;
1110 phydbl *t;
1111 int swap = NO;
1112
1113 for(i=0;i<2*tree->n_otu-1;++i) tree->times->t_rank[i] = i;
1114
1115 t = tree->times->nd_t;
1116
1117 do
1118 {
1119 swap = NO;
1120 for(i=0;i<2*tree->n_otu-2;++i)
1121 {
1122 if(t[tree->times->t_rank[i]] > t[tree->times->t_rank[i+1]]) // Sort in ascending order
1123 {
1124 swap = YES;
1125 buff = tree->times->t_rank[i];
1126 tree->times->t_rank[i] = tree->times->t_rank[i+1];
1127 tree->times->t_rank[i+1] = buff;
1128 }
1129 }
1130 }
1131 while(swap == YES);
1132
1133 }
1134
1135 //////////////////////////////////////////////////////////////
1136 //////////////////////////////////////////////////////////////
1137
1138
TIMES_Set_Calibration(t_tree * tree)1139 void TIMES_Set_Calibration(t_tree *tree)
1140 {
1141 t_cal *cal;
1142 int i;
1143
1144 For(i,2*tree->n_otu-1)
1145 {
1146 tree->times->t_has_prior[i] = NO;
1147 tree->times->t_prior_min[i] = BIG;
1148 tree->times->t_prior_max[i] = BIG;
1149 }
1150
1151 cal = tree->times->a_cal[0];
1152 while(cal)
1153 {
1154 /* if(cal->is_active == YES) */
1155 /* { */
1156 /* tree->times->t_has_prior[cal->node_num] = YES; */
1157 /* tree->times->t_prior_min[cal->node_num] = cal->lower; */
1158 /* tree->times->t_prior_max[cal->node_num] = cal->upper; */
1159 /* } */
1160 cal = cal->next;
1161 }
1162
1163 TIMES_Set_All_Node_Priors(tree);
1164 }
1165
1166 //////////////////////////////////////////////////////////////
1167 //////////////////////////////////////////////////////////////
1168
TIMES_Record_Prior_Times(t_tree * tree)1169 void TIMES_Record_Prior_Times(t_tree *tree)
1170 {
1171 int i;
1172 for(i=0;i<2*tree->n_otu-1;++i)
1173 {
1174 tree->times->t_prior_min_ori[i] = tree->times->t_prior_min[i];
1175 tree->times->t_prior_max_ori[i] = tree->times->t_prior_max[i];
1176 }
1177 }
1178
1179 //////////////////////////////////////////////////////////////
1180 //////////////////////////////////////////////////////////////
1181
TIMES_Reset_Prior_Times(t_tree * tree)1182 void TIMES_Reset_Prior_Times(t_tree *tree)
1183 {
1184 int i;
1185 For(i,2*tree->n_otu-1)
1186 {
1187 tree->times->t_prior_min[i] = tree->times->t_prior_min_ori[i];
1188 tree->times->t_prior_max[i] = tree->times->t_prior_max_ori[i];
1189 }
1190 }
1191
1192 //////////////////////////////////////////////////////////////
1193 //////////////////////////////////////////////////////////////
1194 //////////////////////////////////////////////////////////////
1195 //////////////////////////////////////////////////////////////
1196 // Returns the conditional density of internal node heights
1197 // given the tree height under the Yule model. Uses the order
1198 // statistics 'simplification' as described in Yang and Rannala, 2005.
TIMES_Lk_Yule_Order_Root_Cond(t_tree * tree)1199 phydbl TIMES_Lk_Yule_Order_Root_Cond(t_tree *tree)
1200 {
1201 int j;
1202 phydbl *t,*tf;
1203 t_node *n;
1204 phydbl loglk;
1205 phydbl loglbda;
1206 phydbl lbda;
1207 phydbl *tp_min,*tp_max;
1208 phydbl lower_bound,upper_bound;
1209 phydbl root_height;
1210
1211 tp_min = tree->times->t_prior_min;
1212 tp_max = tree->times->t_prior_max;
1213 tf = tree->times->t_floor;
1214 t = tree->times->nd_t;
1215 n = NULL;
1216 loglbda = log(tree->times->birth_rate);
1217 lbda = tree->times->birth_rate;
1218 lower_bound = -1.;
1219 upper_bound = -1.;
1220 root_height = FABS(tree->times->nd_t[tree->n_root->num]);
1221
1222 /*! Adapted from Equation (6) in T. Stadler's Systematic Biology, 2012 paper with
1223 sampling fraction set to 1 and death rate set to 0. Dropped the 1/(n-1) scaling
1224 factor. */
1225
1226 /* loglk = 0.0; */
1227 /* For(j,2*tree->n_otu-2) */
1228 /* { */
1229 /* n = tree->a_nodes[j]; */
1230 /* lower_bound = MAX(FABS(tf[j]),FABS(tp_max[j])); */
1231 /* upper_bound = MIN(FABS(t[tree->n_root->num]),FABS(tp_min[j])); */
1232
1233 /* if(n->tax == NO) */
1234 /* { */
1235 /* loglk += (loglbda - lbda * FABS(t[j])); */
1236 /* /\* loglk -= log(exp(-lbda*lower_bound) - exp(-lbda*upper_bound)); // incorporate calibration boundaries here. *\/ */
1237 /* } */
1238 /* } */
1239
1240
1241 /*! Adapted from Equation (7) in T. Stadler's Systematic Biology, 2012 paper with
1242 sampling fraction set to 1 and death rate set to 0. */
1243
1244 // Check that each node is within its calibration-derived interval
1245 For(j,2*tree->n_otu-1) if(t[j] < tp_min[j] || t[j] > tp_max[j]) return(-INFINITY);
1246
1247 loglk = 0.0;
1248 For(j,2*tree->n_otu-2)
1249 {
1250 n = tree->a_nodes[j];
1251 lower_bound = MAX(FABS(tf[j]),FABS(tp_max[j]));
1252 upper_bound = MIN(FABS(tp_min[j]),root_height);
1253
1254 if(n->tax == NO)
1255 {
1256 loglk += (loglbda - lbda * FABS(t[j]));
1257 loglk -= log(exp(-lbda*lower_bound) - exp(-lbda*upper_bound)); // incorporate calibration boundaries here.
1258 }
1259
1260 if(isinf(loglk) || isnan(loglk))
1261 {
1262 PhyML_Fprintf(stderr,"\n. Lower bound: %f",lower_bound);
1263 PhyML_Fprintf(stderr,"\n. Upper bound: %f",upper_bound);
1264 PhyML_Fprintf(stderr,"\n. tf: %f tp_max: %f tp_min: %f root: %f",tf[j],tp_max[j],tp_min[j],root_height);
1265 Exit("\n");
1266 }
1267
1268 }
1269
1270 /* lower_bound = MAX(FABS(tf[tree->n_root->num]),FABS(tp_max[tree->n_root->num])); */
1271 /* upper_bound = FABS(tp_min[tree->n_root->num]); */
1272 /* loglk += log(2) + loglbda - 2.*lbda * FABS(t[tree->n_root->num]); */
1273 /* loglk -= log(exp(-2.*lbda*lower_bound) - exp(-2.*lbda*upper_bound)); */
1274
1275
1276 return(loglk);
1277 }
1278
1279 //////////////////////////////////////////////////////////////
1280 //////////////////////////////////////////////////////////////
1281 // Log of prob density of internal node ages conditional on tree height, under
1282 // the birth death process with incomplete sampling.
TIMES_Lk_Birth_Death(int verbose,t_tree * tree)1283 phydbl TIMES_Lk_Birth_Death(int verbose, t_tree *tree)
1284 {
1285 int i,n;
1286 phydbl lnL;
1287 phydbl t,b,d,bmd,logbmd,expmbmd,logb;
1288 phydbl bmin,bmax;
1289 phydbl dmin,dmax;
1290 phydbl lognut1,logp_1t,troot,pt,nut1;
1291
1292 lnL = 0.0;
1293 b = tree->times->birth_rate;
1294 d = tree->times->death_rate;
1295 bmin = tree->times->birth_rate_min;
1296 bmax = tree->times->birth_rate_max;
1297 dmin = tree->times->death_rate_min;
1298 dmax = tree->times->death_rate_min;
1299 bmd = b-d;
1300 logbmd = -1.;
1301 expmbmd = -1.;
1302 logb = -1.;
1303 t = 0.0;
1304 n = tree->n_otu;
1305 troot = fabs(tree->times->nd_t[tree->n_root->num]);
1306 logb = log(b);
1307
1308 if(b < d)
1309 {
1310 tree->times->c_lnL_times = UNLIKELY;
1311 if(verbose) PhyML_Printf("\n. b < d");
1312 return UNLIKELY;
1313 }
1314
1315 // Verify that calibration constraints are satisfied
1316 for(i=0;i<2*tree->n_otu-1;++i)
1317 if(tree->a_nodes[i]->tax == NO)
1318 {
1319 if(tree->times->nd_t[i] < tree->times->t_prior_min[i] ||
1320 tree->times->nd_t[i] > tree->times->t_prior_max[i])
1321 {
1322 tree->times->c_lnL_times = UNLIKELY;
1323 if(verbose)
1324 {
1325 PhyML_Printf("\n. Time outside calibration range : %G [%G,%G]",tree->times->nd_t[i],tree->times->t_prior_min[i],tree->times->t_prior_max[i]);
1326 PhyML_Printf("\n. Clade incriminated: ");
1327 if(tree->a_nodes[i] == tree->n_root)
1328 {
1329 List_Taxa_In_Clade(tree->n_root,tree->n_root->v[1],tree);
1330 List_Taxa_In_Clade(tree->n_root,tree->n_root->v[2],tree);
1331 }
1332 else
1333 {
1334 assert(tree->a_nodes[i]->anc);
1335 List_Taxa_In_Clade(tree->a_nodes[i]->anc,tree->a_nodes[i],tree);
1336 }
1337 PhyML_Printf("\n");
1338 }
1339 return UNLIKELY;
1340 }
1341 }
1342
1343 if(b > bmin && d > dmin && Are_Equal(bmd,0.0,bmin/10.) == NO)
1344 {
1345 logbmd = log(bmd);
1346 expmbmd = exp(d-b);
1347
1348 pt = bmd/(b-d*pow(expmbmd,troot));
1349 nut1 = 1. - pt*pow(expmbmd,troot);
1350 lognut1 = log(nut1);
1351 /* printf("\n. lognut1: %G pt: %G b: %G d: %G exp: %G",nut1,pt,b,d,expmbmd); fflush(NULL); */
1352
1353 for(i=0;i<2*tree->n_otu-1;++i)
1354 if(tree->a_nodes[i]->tax == NO && tree->a_nodes[i] != tree->n_root)
1355 {
1356 t = fabs(tree->times->nd_t[i]);
1357
1358 /* // Equation 3.19 in Tanja Stadler's PhD thesis */
1359 /* lnL += 2*logbmd - bmd*t - 2.*log(b-d*pow(expmbmd,t)); */
1360
1361 // Equation 6 in Yang and Rannala, 1997 with rho=1
1362 logp_1t = 2.*logbmd - 2.*log(b-d*exp((d-b)*t)) + (d-b)*t;
1363 lnL += logb + logp_1t - lognut1;
1364
1365
1366 if(!(lnL > UNLIKELY))
1367 {
1368 PhyML_Printf("\n. logb: %G pt: %G nut1: %G lognut1: %G t: %G logp_1t: %G\n",logb,pt,nut1,lognut1,t,logp_1t);
1369 tree->times->c_lnL_times = UNLIKELY;
1370 return UNLIKELY;
1371 }
1372 }
1373
1374 lnL += LnGamma(n-1);
1375
1376
1377 /* t = FABS(tree->times->nd_t[tree->n_root->num]); */
1378 /* lnL += -bmd*t - log(b-d*pow(expmbmd,t)); */
1379 /* lnL += log(bmd) + (tree->n_otu-1)*log(b) + LnGamma(tree->n_otu+1); */
1380
1381 // Divide joint density p(x(1),...,x(n-1)) by p(x(1)) in order to get
1382 // the conditional density p(x(2),...,x(n-1)|x(1)). The log of the marginal p(x(1))
1383 // (Theorem 3.4.11 in Tanja's PhD thesis) is given below and subtracted to lnL.
1384 /* k = 1; */
1385 /* s = FABS(tree->times->nd_t[tree->n_root->num]); */
1386 /* phydbl p_x1 = (k+1)*Choose(n,k+1)*pow(b,n-k)*pow(bmd,k+2)*exp(-bmd*(k+1)*s)*pow(1.-exp(-bmd*s),n-k-1)/pow(b-d*exp(-bmd*s),n+1); */
1387 /* lnL -= log(p_x1); */
1388 }
1389 else if(b > bmin && d < dmin) // Yule process
1390 {
1391 lognut1 = log(1.-exp(-b*troot));
1392
1393 for(i=0;i<2*tree->n_otu-1;++i)
1394 if(tree->a_nodes[i]->tax == NO && tree->a_nodes[i] != tree->n_root)
1395 {
1396 t = fabs(tree->times->nd_t[i]);
1397 /* // Equation 3.19 in Tanja Stadler's PhD thesis (Yule case) */
1398 /* lnL -= b*t; */
1399
1400 // Equation 6 in Yang and Rannala, 1997 with rho=1
1401 logp_1t = - b*t;
1402 lnL += logb + logp_1t - lognut1;
1403
1404 if(!(lnL > UNLIKELY))
1405 {
1406 PhyML_Printf("\n. lognut1: %G t: %G logp_1t: %G",lognut1,t,logp_1t);
1407 tree->times->c_lnL_times = UNLIKELY;
1408 return UNLIKELY;
1409 }
1410 }
1411
1412 lnL += LnGamma(n-1);
1413
1414 /* t = fabs(tree->times->nd_t[tree->n_root->num]); */
1415 /* lnL -= b*t; */
1416 /* lnL += (tree->n_otu-1)*log(b) + LnGamma(tree->n_otu+1); */
1417
1418
1419 // Divide joint density p(x(1),...,x(n-1)) by p(x(1)) in order to get
1420 // the conditional density p(x(2),...,x(n-1)|x(1)). The log of the marginal p(x(1))
1421 // (Theorem 3.4.11, remark 3.4.12 in Tanja's PhD thesis) is given below and subtracted to lnL.
1422 /* k = 1; */
1423 /* s = FABS(tree->times->nd_t[tree->n_root->num]); */
1424 /* phydbl p_x1 = (k+1)*Choose(n,k+1)*b*pow(exp(b*s)-1.,n-k-1)/exp(b*s*n); */
1425 /* lnL -= log(p_x1); */
1426 }
1427 else if(b < bmin && d > dmin)
1428 {
1429 PhyML_Printf("\n. b: %G bmin: %G d: %G dmin: %G",b,bmin,d,dmin);
1430 tree->times->c_lnL_times = UNLIKELY;
1431 return UNLIKELY;
1432 }
1433 else if(Are_Equal(bmd,0.0,bmin/10.) == YES) // Critical birth-death process
1434 {
1435 logb = log(b);
1436
1437 for(i=0;i<2*tree->n_otu-1;++i)
1438 if(tree->a_nodes[i]->tax == NO && tree->a_nodes[i] != tree->n_root)
1439 {
1440 t = fabs(tree->times->nd_t[i]);
1441
1442 /* // Equation 3.19 in Tanja Stadler's PhD thesis (Critical case) */
1443 /* lnL += logb - 2.*log(1.+b*t); */
1444
1445
1446 // Equation 7 in Yang and Rannala, 1997 with rho=1
1447 lnL += log((1.+d)/pow(1.+d*t,2));
1448
1449 if(!(lnL > UNLIKELY))
1450 {
1451 PhyML_Printf("\n. logb: %G t: %G",logb,t);
1452 tree->times->c_lnL_times = UNLIKELY;
1453 return UNLIKELY;
1454 }
1455 }
1456
1457 lnL += LnGamma(n-1);
1458
1459 /* t = fabs(tree->times->nd_t[tree->n_root->num]); */
1460 /* lnL -= log(b*t); */
1461 /* lnL += LnGamma(tree->n_otu+1); */
1462
1463 // Divide joint density p(x(1),...,x(n-1)) by p(x(1)) in order to get
1464 // the conditional density p(x(2),...,x(n-1)|x(1)). The log of the marginal p(x(1))
1465 // (Theorem 3.4.11, remark 3.4.12 in Tanja's PhD thesis) is given below and subtracted to lnL.
1466 /* k = 1; */
1467 /* s = fabs(tree->times->nd_t[tree->n_root->num]); */
1468 /* phydbl p_x1 = (k+1)*Choose(n,k+1)*pow(b,n-k)*pow(s,n-k-1)/pow(1.+b*s,n+1); */
1469 /* lnL -= log(p_x1); */
1470 }
1471 else if(b < bmin && d < dmin) // Birth and death rates are below their limits
1472 {
1473 PhyML_Fprintf(stderr,"\n. b: %G bmin: %G d: %G dmin: %G",b,bmin,d,dmin);
1474 tree->times->c_lnL_times = UNLIKELY;
1475 return -INFINITY;
1476 }
1477 else if(b > bmax && d > dmax)
1478 {
1479 PhyML_Fprintf(stderr,"\n. b: %G bmax: %G d: %G dmax: %G",b,bmax,d,dmax);
1480 tree->times->c_lnL_times = UNLIKELY;
1481 return -INFINITY;
1482 }
1483 else
1484 {
1485 assert(FALSE);
1486 }
1487
1488 if(isnan(lnL) || isinf(fabs(lnL)) || !(lnL > UNLIKELY))
1489 {
1490 PhyML_Fprintf(stderr,"\n. lnL times: %f",lnL);
1491 tree->times->c_lnL_times = UNLIKELY;
1492 return UNLIKELY;
1493 }
1494
1495 tree->times->c_lnL_times = lnL;
1496
1497 return(lnL);
1498 }
1499
1500 //////////////////////////////////////////////////////////////
1501 //////////////////////////////////////////////////////////////
1502
1503 // Generate a subtree including all taxa in tax_list. The age of the root of that
1504 // subtree is t_mrca. All nodes in the subtree are thus younger than that.
TIMES_Connect_List_Of_Taxa(t_node ** tax_list,int list_size,phydbl t_mrca,phydbl * times,int * nd_num,t_tree * mixt_tree)1505 void TIMES_Connect_List_Of_Taxa(t_node **tax_list, int list_size, phydbl t_mrca, phydbl *times, int *nd_num, t_tree *mixt_tree)
1506 {
1507 phydbl t_upper_bound, t_lower_bound,up,lo;
1508 int i,j,n_anc,*permut,rand_idx;
1509 t_node *n,**anc,*new_mrca;
1510
1511 // Calibration of a tip node
1512 if(list_size == 1)
1513 {
1514 times[tax_list[0]->num] = t_mrca;
1515 return;
1516 }
1517 else
1518 {
1519 t_lower_bound = t_mrca;
1520 t_upper_bound = +INFINITY;
1521 n = NULL;
1522 anc = NULL;
1523 new_mrca = NULL;
1524 permut = NULL;
1525
1526 // Find the upper bound for all the new node ages that
1527 // will be created in this function
1528 for(i=0;i<list_size;i++)
1529 {
1530 n = tax_list[i];
1531 while(n->v[0] != NULL) n = n->v[0];
1532 if(times[n->num] < t_upper_bound) t_upper_bound = times[n->num];
1533 }
1534
1535 if(t_upper_bound < t_lower_bound)
1536 {
1537 PhyML_Printf("\n. upper: %f lower: %f t_mrca: %f\n",t_upper_bound,t_lower_bound,t_mrca);
1538 assert(FALSE);
1539 }
1540
1541 // Get the list of current mrcas to all taxa in tax_list. There should be
1542 // at least one of these
1543 n_anc = 0;
1544 for(i=0;i<list_size;i++)
1545 {
1546 n = tax_list[i];
1547 while(n->v[0] != NULL) n = n->v[0];
1548 for(j=0;j<n_anc;j++) if(anc[j] == n) break;
1549 if(j == n_anc)
1550 {
1551 if(n_anc == 0) anc = (t_node **)mCalloc(1,sizeof(t_node *));
1552 else anc = (t_node **)mRealloc(anc,n_anc+1,sizeof(t_node *));
1553 anc[n_anc] = n;
1554 n_anc++;
1555 }
1556 }
1557
1558 /* printf("\n. n_anc: %d",n_anc); */
1559 /* for(i=0;i<n_anc;i++) PhyML_Printf("\n. anc: %d",anc[i]->num); */
1560
1561 if(n_anc == 1) // All the nodes in tax_list are already connected. Bail out.
1562 {
1563 Free(anc);
1564 return;
1565 }
1566
1567 // Connect randomly and set ages
1568 permut = Permutate(n_anc);
1569 /* permut = (int *)mCalloc(n_anc,sizeof(int)); */
1570 /* for(i=0;i<n_anc;i++) permut[i] = i; */
1571 i = 0;
1572 do
1573 {
1574 new_mrca = mixt_tree->a_nodes[*nd_num];
1575 anc[permut[i]]->v[0] = new_mrca;
1576 anc[permut[i+1]]->v[0] = new_mrca;
1577 new_mrca->v[1] = anc[permut[i]];
1578 new_mrca->v[2] = anc[permut[i+1]];
1579 new_mrca->v[0] = NULL;
1580 up = MIN(times[new_mrca->v[1]->num],times[new_mrca->v[2]->num]);
1581 lo = up - (up - t_mrca)/5.;
1582 times[new_mrca->num] = Uni()*(up - lo) + lo;
1583
1584 /* times[new_mrca->num] = t_upper_bound - ((phydbl)(i+1.)/n_anc)*(t_upper_bound - t_lower_bound); */
1585
1586 /* printf("\n. new_mrca->num: %d time: %f [%f %f] t_mrca: %f %d connect to %d %d %d [%f %f]", */
1587 /* new_mrca->num, */
1588 /* times[new_mrca->num], */
1589 /* t_lower_bound, */
1590 /* t_upper_bound, */
1591 /* t_mrca, */
1592 /* new_mrca->num, */
1593 /* new_mrca->v[0] ? new_mrca->v[0]->num : -1, */
1594 /* new_mrca->v[1] ? new_mrca->v[1]->num : -1, */
1595 /* new_mrca->v[2] ? new_mrca->v[2]->num : -1, */
1596 /* times[new_mrca->v[1]->num], */
1597 /* times[new_mrca->v[2]->num] */
1598 /* ); */
1599 /* fflush(NULL); */
1600
1601 anc[permut[i+1]] = new_mrca;
1602
1603 rand_idx = Rand_Int(i+1,n_anc-1);
1604 n = anc[permut[i+1]];
1605 anc[permut[i+1]] = anc[permut[rand_idx]];
1606 anc[permut[rand_idx]] = n;
1607
1608
1609 i++;
1610 (*nd_num) += 1;
1611 if(n_anc == i+1) { times[new_mrca->num] = t_mrca; break; }
1612 }
1613 while(1);
1614
1615 Free(permut);
1616 Free(anc);
1617 }
1618 }
1619 //////////////////////////////////////////////////////////////
1620 //////////////////////////////////////////////////////////////
1621 // Generate a random rooted tree with node ages fullfiling the
1622 // time constraints defined in the list of calibration cal_list
TIMES_Randomize_Tree_With_Time_Constraints(t_cal * cal_list,t_tree * mixt_tree)1623 void TIMES_Randomize_Tree_With_Time_Constraints(t_cal *cal_list, t_tree *mixt_tree)
1624 {
1625 t_node **tips,**nd_list;
1626 phydbl *times,*cal_times,time_oldest_cal;
1627 int i,j,nd_num,*cal_ordering,n_cal,swap,list_size,tmp,orig_is_mixt_tree,repeat,n_max_repeats,tip_num,*no_cal_tip_num,n_no_cal_tip_num,*permut,*tip_is_in_cal_clad;
1628 t_cal *cal;
1629 t_clad *clade;
1630
1631 assert(mixt_tree->rates);
1632
1633 if(mixt_tree->mod->s_opt->opt_topo == NO)
1634 {
1635 PhyML_Fprintf(stderr,"\n. Fixing the tree topology is only allowed when calibrating tip nodes only.");
1636 PhyML_Fprintf(stderr,"\n. You are most likely calibrating here the MRCA of at least one clade with more than two tips.");
1637 PhyML_Fprintf(stderr,"\n. It is difficult to set the age of that clade within the limit of the calibration constraints,");
1638 PhyML_Fprintf(stderr,"\n. and fix the tree topology at the same time. Please contact me (guindon@lirmm.fr) for a more");
1639 PhyML_Fprintf(stderr,"\n. detailed diagnostic.");
1640 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
1641 }
1642
1643 tips = (t_node **)mCalloc(mixt_tree->n_otu,sizeof(t_node *));
1644 nd_list = (t_node **)mCalloc(2*mixt_tree->n_otu-1,sizeof(t_node *));
1645 tip_is_in_cal_clad = (int *)mCalloc(mixt_tree->n_otu,sizeof(t_node *));
1646
1647 times = mixt_tree->times->nd_t;
1648 orig_is_mixt_tree = mixt_tree->is_mixt_tree;
1649 mixt_tree->is_mixt_tree = NO;
1650 n_max_repeats = 1000;
1651 cal = NULL;
1652 clade = NULL;
1653
1654
1655 // Is tip in calibrated clade or not?
1656 for(i=0;i<mixt_tree->n_otu;++i)
1657 {
1658 cal = cal_list;
1659 clade = NULL;
1660 while(cal != NULL)
1661 {
1662 while(cal && cal->clade_list == NULL) cal = cal->next;
1663 if(cal == NULL) break;
1664 clade = cal->clade_list[cal->current_clade_idx];
1665 for(j=0;j<clade->n_tax;++j) if(clade->tip_list[j] == mixt_tree->a_nodes[i]) break;
1666 if(j != clade->n_tax) break;
1667 cal = cal->next;
1668 }
1669
1670 if(cal != NULL) tip_is_in_cal_clad[i] = YES;
1671 else tip_is_in_cal_clad[i] = NO;
1672 }
1673
1674 // List node indices that are not in any calibration set
1675 no_cal_tip_num = NULL;
1676 n_no_cal_tip_num = 0;
1677 for(i=0;i<mixt_tree->n_otu;++i)
1678 {
1679 cal = cal_list;
1680 while(cal != NULL)
1681 {
1682 while(cal && cal->clade_list == NULL) cal = cal->next;
1683 if(cal == NULL) break;
1684 clade = cal->clade_list[cal->current_clade_idx];
1685 for(j=0;j<clade->n_tax;++j) if(clade->tip_list[j] == mixt_tree->a_nodes[i]) break;
1686 if(j != clade->n_tax) break;
1687 cal = cal->next;
1688 }
1689
1690 if(cal == NULL)
1691 {
1692 if(n_no_cal_tip_num == 0) no_cal_tip_num = (int *)mCalloc(1,sizeof(int));
1693 else no_cal_tip_num = (int *)mRealloc(no_cal_tip_num,n_no_cal_tip_num+1,sizeof(int));
1694 no_cal_tip_num[n_no_cal_tip_num] = i;
1695 n_no_cal_tip_num++;
1696 }
1697 }
1698
1699
1700 for(repeat=0;repeat<n_max_repeats;++repeat)
1701 {
1702 nd_num = mixt_tree->n_otu;
1703
1704 /* PhyML_Printf("\n\n. Repeat %d",repeat); */
1705
1706 for(i=0;i<mixt_tree->n_otu;++i) tips[i] = mixt_tree->a_nodes[i];
1707 for(i=0;i<mixt_tree->n_otu;++i) mixt_tree->a_nodes[i]->v[0] = NULL;
1708
1709
1710 // Set a time for each calibration
1711 cal_times = (phydbl *)mCalloc(1,sizeof(phydbl));
1712 time_oldest_cal = 0.0;
1713 n_cal = 0;
1714 cal = cal_list;
1715 while(cal != NULL)
1716 {
1717 if(cal->is_primary == YES)
1718 {
1719 if(n_cal > 0) cal_times = (phydbl *)mRealloc(cal_times,n_cal+1,sizeof(phydbl));
1720 cal_times[n_cal] = Uni()*(cal->upper - cal->lower) + cal->lower;
1721 if(cal_times[n_cal] < time_oldest_cal) time_oldest_cal = cal_times[n_cal];
1722 n_cal++;
1723 }
1724 cal = cal->next;
1725 }
1726
1727 cal_ordering = (int *)mCalloc(n_cal,sizeof(int));
1728 for(i=0;i<n_cal;++i) cal_ordering[i] = i;
1729
1730 // Sort calibration times from youngest to oldest
1731 do
1732 {
1733 swap = NO;
1734 for(i=0;i<n_cal-1;i++)
1735 {
1736 if(cal_times[cal_ordering[i]] < cal_times[cal_ordering[i+1]])
1737 {
1738 tmp = cal_ordering[i+1];
1739 cal_ordering[i+1] = cal_ordering[i];
1740 cal_ordering[i] = tmp;
1741 swap = YES;
1742 }
1743 }
1744 }
1745 while(swap == YES);
1746
1747 for(i=0;i<n_cal-1;i++) assert(!(cal_times[cal_ordering[i]] < cal_times[cal_ordering[i+1]]));
1748
1749
1750 // Connect taxa that appear in every primary calibrations
1751 for(i=0;i<n_cal;++i)
1752 {
1753 cal = cal_list;
1754 j = 0;
1755 while(j != cal_ordering[i])
1756 {
1757 if(cal->is_primary == YES) j++;
1758 cal = cal->next;
1759 assert(cal);
1760 }
1761
1762 list_size = 0;
1763
1764 /* printf("\n. n_cal: %d n_no_cal_tip_num: %d [%d %d] n_otu: %d", */
1765 /* n_cal, */
1766 /* n_no_cal_tip_num, */
1767 /* n_no_cal_tip_num-1, */
1768 /* (int)(2*n_no_cal_tip_num/(n_cal)), */
1769 /* mixt_tree->n_otu); fflush(NULL); */
1770
1771 // Add some taxa that are not in any calibration set
1772 if(n_no_cal_tip_num > 0)
1773 {
1774 permut = Permutate(n_no_cal_tip_num);
1775 j = 0;
1776 do
1777 {
1778 tip_num = no_cal_tip_num[permut[j]];
1779 if(tips[tip_num]->v[0] == NULL)
1780 {
1781 nd_list[list_size] = tips[tip_num];
1782 /* PhyML_Printf("\n# %s",tips[tip_num]->name); */
1783 list_size++;
1784 assert(list_size <= mixt_tree->n_otu);
1785 }
1786 }
1787 while(j++ < MIN(n_no_cal_tip_num-1,(int)(2*n_no_cal_tip_num/(n_cal))));
1788 Free(permut);
1789 }
1790
1791
1792 // Add all the taxa that are in the calibration set of cal
1793 // This should be done here so that the last node in nd_list
1794 // belongs to the taxa in the calibration
1795
1796 if(cal->clade_list != NULL)
1797 {
1798 clade = cal->clade_list[cal->current_clade_idx];
1799 for(j=0;j<clade->n_tax;j++)
1800 {
1801 nd_list[list_size] = tips[clade->tip_list[j]->num];
1802 list_size++;
1803 assert(list_size <= mixt_tree->n_otu);
1804 }
1805
1806 TIMES_Connect_List_Of_Taxa(nd_list,
1807 list_size,
1808 cal_times[cal_ordering[i]],
1809 times,
1810 &nd_num,
1811 mixt_tree);
1812 }
1813 }
1814
1815 Free(cal_times);
1816 Free(cal_ordering);
1817
1818
1819 // Connect all remaining taxa
1820 for(i=0;i<mixt_tree->n_otu;i++) nd_list[i] = NULL;
1821 list_size = 0;
1822 for(i=0;i<mixt_tree->n_otu;++i)
1823 {
1824 if(tip_is_in_cal_clad[tips[i]->num] == NO) // Tip is not in a calibrated clade
1825 {
1826 nd_list[list_size] = tips[i];
1827 list_size++;
1828 assert(list_size <= mixt_tree->n_otu);
1829 }
1830 }
1831
1832 cal = cal_list;
1833 do
1834 {
1835 while(cal && cal->clade_list == NULL) cal = cal->next;
1836 if(cal == NULL) break;
1837 clade = cal->clade_list[cal->current_clade_idx];
1838 nd_list[list_size] = clade->tip_list[0];
1839 list_size++;
1840 assert(list_size <= 2*mixt_tree->n_otu-1);
1841 cal = cal->next;
1842 }
1843 while(cal);
1844
1845 /* for(i=0;i<list_size;i++) printf("\n# To connect: %d lower: %f",nd_list[i]->num,time_oldest_cal); */
1846
1847
1848 TIMES_Connect_List_Of_Taxa(nd_list,
1849 list_size,
1850 10.*time_oldest_cal-1.E-3, // 1.E-3 required in case time_oldest_cal = 0.
1851 times,
1852 &nd_num,
1853 mixt_tree);
1854
1855 // Adding root node
1856 mixt_tree->n_root = mixt_tree->a_nodes[2*mixt_tree->n_otu-2];
1857 mixt_tree->n_root->v[1]->v[0] = mixt_tree->n_root->v[2];
1858 mixt_tree->n_root->v[2]->v[0] = mixt_tree->n_root->v[1];
1859 mixt_tree->n_root->anc = NULL;
1860
1861 Connect_Edges_To_Nodes_Serial(mixt_tree);
1862
1863
1864 // Adding root edge
1865 for(i=0;i<2*mixt_tree->n_otu-3;++i)
1866 {
1867 if(((mixt_tree->a_edges[i]->left == mixt_tree->n_root->v[1]) || (mixt_tree->a_edges[i]->rght == mixt_tree->n_root->v[1])) &&
1868 ((mixt_tree->a_edges[i]->left == mixt_tree->n_root->v[2]) || (mixt_tree->a_edges[i]->rght == mixt_tree->n_root->v[2])))
1869 {
1870 mixt_tree->e_root = mixt_tree->a_edges[i];
1871 break;
1872 }
1873 }
1874
1875 Update_Ancestors(mixt_tree->n_root,mixt_tree->n_root->v[2],mixt_tree);
1876 Update_Ancestors(mixt_tree->n_root,mixt_tree->n_root->v[1],mixt_tree);
1877 DATE_Assign_Primary_Calibration(mixt_tree);
1878 DATE_Update_T_Prior_MinMax(mixt_tree);
1879
1880
1881 /* for(int i=0;i<mixt_tree->times->n_cal;i++) */
1882 /* { */
1883 /* int idx; */
1884 /* cal = mixt_tree->times->a_cal[i]; */
1885
1886 /* if(cal->clade_list != NULL) */
1887 /* { */
1888 /* clade = cal->clade_list[cal->current_clade_idx]; */
1889 /* idx = Find_Clade(clade->tax_list,clade->n_tax,mixt_tree); */
1890 /* PhyML_Printf("\n. Calibration %3d | Node number to which calibration applies is [%d]",i,idx); */
1891 /* PhyML_Printf("\n. Calibration %3d | Lower bound set to: %15f time units.",i,mixt_tree->times->a_cal[i]->lower); */
1892 /* PhyML_Printf("\n. Calibration %3d | Upper bound set to: %15f time units.",i,mixt_tree->times->a_cal[i]->upper); */
1893 /* PhyML_Printf("\n. Calibration %3d | t_prior_min: %G t_prior_max: %G",i,mixt_tree->times->t_prior_min[idx],mixt_tree->times->t_prior_max[idx]); */
1894 /* PhyML_Printf("\n. Calibration %3d | Time set to %G",i,mixt_tree->times->nd_t[idx]); */
1895 /* } */
1896 /* } */
1897
1898 if(!DATE_Check_Calibration_Constraints(mixt_tree))
1899 {
1900 /* PhyML_Printf("\n. Could not generate tree (DATE_Check_Calibration_Constraints)\n\n"); */
1901 }
1902 else if(!DATE_Check_Time_Constraints(mixt_tree))
1903 {
1904 /* PhyML_Printf("\n. Could not generate tree (DATE_Check_Time_Constraints)\n\n"); */
1905 }
1906 else break; // Tree successfully generated
1907 }
1908
1909
1910 if(repeat == n_max_repeats)
1911 {
1912 PhyML_Fprintf(stderr,"\n\n");
1913 PhyML_Fprintf(stderr,"\n. A random tree satisfying the calibration constraints provided");
1914 PhyML_Fprintf(stderr,"\n. could not be generated. It probably means that there are some");
1915 PhyML_Fprintf(stderr,"\n. inconsistencies in the calibration data. For instance, the calibration");
1916 PhyML_Fprintf(stderr,"\n. time interval for the MRCA of a clade with taxa {X,Y} (noted as [a,b])");
1917 PhyML_Fprintf(stderr,"\n. cannot be strictly older than the interval corresponding to taxa ");
1918 PhyML_Fprintf(stderr,"\n. {X,Z,Y} (noted as [c,d]), i.e., b cannot be smaller (older) than c. ");
1919 PhyML_Fprintf(stderr,"\n. Also, please remember that the present time corresponds to a time");
1920 PhyML_Fprintf(stderr,"\n. value equal to zero and past events have negative time values.");
1921
1922 PhyML_Fprintf(stderr,"\n\n");
1923 if(!DATE_Check_Calibration_Constraints(mixt_tree))
1924 {
1925 PhyML_Fprintf(stderr,"\n. Could not generate tree (DATE_Check_Calibration_Constraints)\n\n");
1926 }
1927 else if(!DATE_Check_Time_Constraints(mixt_tree))
1928 {
1929 PhyML_Fprintf(stderr,"\n. Could not generate tree (DATE_Check_Time_Constraints)\n\n");
1930 }
1931
1932 Exit("\n");
1933 }
1934
1935 /* assert(i != 2*mixt_tree->n_otu-3); */
1936
1937 mixt_tree->is_mixt_tree = orig_is_mixt_tree;
1938
1939 Free(tips);
1940 Free(nd_list);
1941 Free(no_cal_tip_num);
1942 Free(tip_is_in_cal_clad);
1943 }
1944
1945 //////////////////////////////////////////////////////////////
1946 //////////////////////////////////////////////////////////////
1947
TIMES_Check_Node_Height_Ordering(t_tree * tree)1948 int TIMES_Check_Node_Height_Ordering(t_tree *tree)
1949 {
1950 if(!TIMES_Check_Node_Height_Ordering_Post(tree->n_root,tree->n_root->v[1],tree)) return NO;
1951 if(!TIMES_Check_Node_Height_Ordering_Post(tree->n_root,tree->n_root->v[2],tree)) return NO;
1952 return YES;
1953 }
1954
1955 //////////////////////////////////////////////////////////////
1956 //////////////////////////////////////////////////////////////
1957
TIMES_Check_Node_Height_Ordering_Post(t_node * a,t_node * d,t_tree * tree)1958 int TIMES_Check_Node_Height_Ordering_Post(t_node *a, t_node *d, t_tree *tree)
1959 {
1960
1961 if(d->anc != a)
1962 {
1963 PhyML_Printf("\n. d=%d d->anc=%d a=%d root=%d",d->num,d->anc->num,a->num,tree->n_root->num);
1964 return NO;
1965 }
1966 if(tree->times->nd_t[d->num] < tree->times->nd_t[a->num])
1967 {
1968 PhyML_Printf("\n. a->t = %f [num:%d] d->t %f [num:%d]",
1969 tree->times->nd_t[a->num],
1970 a->num,
1971 tree->times->nd_t[d->num],
1972 d->num);
1973 return NO;
1974 }
1975 if(d->tax == YES) return YES;
1976 else
1977 {
1978 int i;
1979 for(i=0;i<3;i++)
1980 {
1981 if(d->v[i] != a && d->b[i] != tree->e_root)
1982 if(!TIMES_Check_Node_Height_Ordering_Post(d,d->v[i],tree))
1983 return NO;
1984 }
1985 }
1986 return YES;
1987 }
1988
1989 //////////////////////////////////////////////////////////////
1990 //////////////////////////////////////////////////////////////
1991
1992 /*
1993 2n-2 edge lengths (in a rooted tree), n-1 internal node
1994 ages. Return \sum_i (l_i - delta t_i)^2 where delta t_i is the
1995 calendar time elapsed along the edge of length l_i
1996 */
TIMES_Least_Square_Criterion(t_tree * tree)1997 phydbl TIMES_Least_Square_Criterion(t_tree *tree)
1998 {
1999 phydbl score;
2000 score = 0.0;
2001 assert(tree->n_root);
2002 assert(tree->rates);
2003 TIMES_Pre_Least_Square_Criterion(tree->n_root,tree->n_root->v[1],tree->n_root->b[1],&score,tree);
2004 TIMES_Pre_Least_Square_Criterion(tree->n_root,tree->n_root->v[2],tree->n_root->b[2],&score,tree);
2005 return(score);
2006 }
2007
2008 //////////////////////////////////////////////////////////////
2009 //////////////////////////////////////////////////////////////
2010
TIMES_Pre_Least_Square_Criterion(t_node * a,t_node * d,t_edge * b,phydbl * score,t_tree * tree)2011 void TIMES_Pre_Least_Square_Criterion(t_node *a, t_node *d, t_edge *b, phydbl *score, t_tree *tree)
2012 {
2013 int i;
2014
2015 (*score) += pow(b->l->v - fabs(tree->times->nd_t[a->num] + tree->times->nd_t[d->num])*tree->rates->clock_r,2);
2016
2017 if(d->tax) return;
2018
2019 for(i=0;i<3;++i)
2020 if(d->v[i] != a && d->b[i] != tree->e_root)
2021 TIMES_Pre_Least_Square_Criterion(d,d->v[i],d->b[i],score,tree);
2022 }
2023
2024 //////////////////////////////////////////////////////////////
2025 //////////////////////////////////////////////////////////////
2026
TIMES_Randomize_Node_Ages(t_tree * tree)2027 void TIMES_Randomize_Node_Ages(t_tree *tree)
2028 {
2029 assert(tree->n_root);
2030 assert(tree->rates);
2031 TIMES_Post_Randomize_Node_Ages(tree->n_root,tree->n_root->v[1],tree);
2032 TIMES_Post_Randomize_Node_Ages(tree->n_root,tree->n_root->v[2],tree);
2033 tree->times->nd_t[tree->n_root->num] =
2034 MIN(tree->times->nd_t[tree->n_root->v[1]->num],
2035 tree->times->nd_t[tree->n_root->v[2]->num])
2036 - Rexp(1.);
2037 /* PhyML_Printf("\n. RAND TIMES node %3d %15f",tree->n_root->num,tree->times->nd_t[tree->n_root->num]); */
2038 }
2039
2040 //////////////////////////////////////////////////////////////
2041 //////////////////////////////////////////////////////////////
2042
TIMES_Post_Randomize_Node_Ages(t_node * a,t_node * d,t_tree * tree)2043 void TIMES_Post_Randomize_Node_Ages(t_node *a, t_node *d, t_tree *tree)
2044 {
2045 if(d->tax == YES) return;
2046 else
2047 {
2048 int i,dir1,dir2;
2049
2050 for(i=0;i<3;++i)
2051 if(d->v[i] != a && d->b[i] != tree->e_root)
2052 TIMES_Post_Randomize_Node_Ages(d,d->v[i],tree);
2053
2054 dir1 = dir2 = -1;
2055 for(i=0;i<3;++i)
2056 if(d->v[i] != a && d->b[i] != tree->e_root)
2057 {
2058 if(dir1 < 0) dir1 = i;
2059 else dir2 = i;
2060 }
2061
2062 tree->times->nd_t[d->num] =
2063 MIN(tree->times->nd_t[d->v[dir1]->num],
2064 tree->times->nd_t[d->v[dir2]->num])
2065 - Rexp(1.);
2066 /* PhyML_Printf("\n. RAND TIMES node %3d %15f",d->num,tree->times->nd_t[d->num]); */
2067 }
2068 }
2069
2070
2071 //////////////////////////////////////////////////////////////
2072 //////////////////////////////////////////////////////////////
2073
TIMES_Calibrations_Apply_To_Tips_Only(t_tree * tree)2074 int TIMES_Calibrations_Apply_To_Tips_Only(t_tree *tree)
2075 {
2076 t_cal *cal;
2077 t_clad *clade;
2078
2079 cal = tree->times->a_cal[0];
2080 assert(cal);
2081 clade = NULL;
2082
2083 do
2084 {
2085 while(cal && cal->clade_list == NULL) cal = cal->next;
2086 if(cal == NULL) break;
2087 clade = cal->clade_list[cal->current_clade_idx];
2088 if(clade && clade->n_tax > 1) break;
2089 cal = cal->next;
2090 }
2091 while(cal);
2092
2093 if(cal == NULL) return YES;
2094
2095 return NO;
2096 }
2097
2098 //////////////////////////////////////////////////////////////
2099 //////////////////////////////////////////////////////////////
2100
TIMES_Randomize_Tip_Times_Given_Calibrations(t_tree * tree)2101 void TIMES_Randomize_Tip_Times_Given_Calibrations(t_tree *tree)
2102 {
2103
2104 t_cal *cal;
2105 t_clad *clade;
2106
2107 cal = tree->times->a_cal[0];
2108 assert(cal);
2109
2110 do
2111 {
2112 clade = cal->clade_list[cal->current_clade_idx];
2113
2114 if(clade->n_tax == 1)
2115 {
2116 assert(clade->target_nd->tax == YES);
2117 tree->times->nd_t[clade->target_nd->num] = Uni()*(cal->upper - cal->lower) + cal->lower;
2118 }
2119
2120 cal = cal->next;
2121 }
2122 while(cal);
2123 }
2124
2125 //////////////////////////////////////////////////////////////
2126 //////////////////////////////////////////////////////////////
2127
TIMES_Time_To_Bl(t_tree * tree)2128 void TIMES_Time_To_Bl(t_tree *tree)
2129 {
2130 TIMES_Time_To_Bl_Pre(tree->n_root,tree->n_root->v[1],tree->n_root->b[1],tree);
2131 TIMES_Time_To_Bl_Pre(tree->n_root,tree->n_root->v[2],tree->n_root->b[2],tree);
2132 tree->n_root->b[1]->l->v = tree->times->nd_t[tree->n_root->v[1]->num] - tree->times->nd_t[tree->n_root->num];
2133 tree->n_root->b[2]->l->v = tree->times->nd_t[tree->n_root->v[2]->num] - tree->times->nd_t[tree->n_root->num];
2134 tree->e_root->l->v = tree->n_root->b[1]->l->v + tree->n_root->b[2]->l->v;
2135 }
2136
2137 //////////////////////////////////////////////////////////////
2138 //////////////////////////////////////////////////////////////
2139
TIMES_Time_To_Bl_Pre(t_node * a,t_node * d,t_edge * b,t_tree * tree)2140 void TIMES_Time_To_Bl_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
2141 {
2142 int i;
2143
2144 b->l->v = tree->times->nd_t[d->num] - tree->times->nd_t[a->num];
2145
2146 if(d->tax) return;
2147 else
2148 {
2149 for(i=0;i<3;i++)
2150 if((d->v[i] != a) && (d->b[i] != tree->e_root))
2151 TIMES_Time_To_Bl_Pre(d,d->v[i],d->b[i],tree);
2152 }
2153 }
2154
2155 //////////////////////////////////////////////////////////////
2156 //////////////////////////////////////////////////////////////
2157
TIMES_Tree_Length(t_tree * tree)2158 phydbl TIMES_Tree_Length(t_tree *tree)
2159 {
2160 phydbl sum;
2161
2162 assert(tree->rates);
2163 assert(tree->e_root);
2164
2165 sum = 0.0;
2166 for(int i=0;i<2*tree->n_otu-1;++i)
2167 if(tree->a_edges[i] != tree->e_root)
2168 sum += fabs(tree->times->nd_t[tree->a_edges[i]->left->num]-
2169 tree->times->nd_t[tree->a_edges[i]->rght->num]);
2170
2171 sum += fabs(tree->times->nd_t[tree->n_root->num] - tree->times->nd_t[tree->n_root->v[1]->num]);
2172 sum += fabs(tree->times->nd_t[tree->n_root->num] - tree->times->nd_t[tree->n_root->v[2]->num]);
2173
2174 return(sum);
2175 }
2176
2177 //////////////////////////////////////////////////////////////
2178 //////////////////////////////////////////////////////////////
2179
TIMES_Bl_To_Times(t_tree * tree)2180 void TIMES_Bl_To_Times(t_tree *tree)
2181 {
2182 t_node *v1,*v2;
2183 int dir1,dir2;
2184 phydbl t1,t2;
2185
2186 TIMES_Bl_To_Times_Post(tree->n_root,tree->n_root->v[1],tree->n_root->b[1],tree);
2187 TIMES_Bl_To_Times_Post(tree->n_root,tree->n_root->v[2],tree->n_root->b[2],tree);
2188
2189 dir1 = 1;
2190 dir2 = 2;
2191
2192 v1 = tree->n_root->v[dir1];
2193 v2 = tree->n_root->v[dir2];
2194
2195 t1 = tree->times->nd_t[v1->num] - MIXT_Get_Mean_Edge_Len(tree->n_root->b[1],tree) / (tree->rates->clock_r * tree->rates->br_r[v1->num]);
2196 t2 = tree->times->nd_t[v2->num] - MIXT_Get_Mean_Edge_Len(tree->n_root->b[2],tree) / (tree->rates->clock_r * tree->rates->br_r[v2->num]);
2197
2198 if(Are_Equal(t1,t2,1.E-6) == NO)
2199 {
2200 PhyML_Fprintf(stderr,"\n. It looks as if the edge lengths suplied do not define an ultrametric tree.");
2201 PhyML_Fprintf(stderr,"\n. Please amend these lengths so as it becomes straightforward to transform your tree");
2202 PhyML_Fprintf(stderr,"\n. into a time-tree. Please contact me (guindon@lirmm.fr) for more information.");
2203 PhyML_Fprintf(stderr,"\n. l1: %f l2: %f",MIXT_Get_Mean_Edge_Len(tree->n_root->b[1],tree),MIXT_Get_Mean_Edge_Len(tree->n_root->b[2],tree));
2204 PhyML_Fprintf(stderr,"\n. t1: %f t2: %f",tree->times->nd_t[v1->num],tree->times->nd_t[v2->num]);
2205 PhyML_Fprintf(stderr,"\n. rr1: %f rr2: %f",tree->rates->br_r[v1->num],tree->rates->br_r[v2->num]);
2206 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
2207 }
2208
2209 tree->times->nd_t[tree->n_root->num] = t1;
2210
2211
2212 }
2213
2214 //////////////////////////////////////////////////////////////
2215 //////////////////////////////////////////////////////////////
2216
TIMES_Bl_To_Times_Post(t_node * a,t_node * d,t_edge * b,t_tree * tree)2217 void TIMES_Bl_To_Times_Post(t_node *a, t_node *d, t_edge *b, t_tree *tree)
2218 {
2219 if(d->tax == YES) return;
2220 else
2221 {
2222 t_node *v1,*v2;
2223 int dir1,dir2;
2224 phydbl t1,t2;
2225
2226 dir1 = dir2 = -1;
2227 for(int i=0;i<3;i++)
2228 if((d->v[i] != a) && (d->b[i] != tree->e_root))
2229 {
2230 TIMES_Bl_To_Times_Post(d,d->v[i],d->b[i],tree);
2231 if(dir1 < 0) dir1 = i;
2232 else dir2 = i;
2233 }
2234
2235 v1 = d->v[dir1];
2236 v2 = d->v[dir2];
2237
2238 t1 = tree->times->nd_t[v1->num] - MIXT_Get_Mean_Edge_Len(d->b[dir1],tree) / (tree->rates->clock_r * tree->rates->br_r[v1->num]);
2239 t2 = tree->times->nd_t[v2->num] - MIXT_Get_Mean_Edge_Len(d->b[dir2],tree) / (tree->rates->clock_r * tree->rates->br_r[v2->num]);
2240
2241 if(Are_Equal(t1,t2,1.E-6) == NO)
2242 {
2243 PhyML_Fprintf(stderr,"\n. It looks at if the edge lengths suplied do not define an ultrametric tree.");
2244 PhyML_Fprintf(stderr,"\n. Please amend these lengths so as it becomes straightforward to transform your tree");
2245 PhyML_Fprintf(stderr,"\n. into a time-tree.");
2246 PhyML_Fprintf(stderr,"\n. l1: %f l2: %f",MIXT_Get_Mean_Edge_Len(d->b[dir1],tree),MIXT_Get_Mean_Edge_Len(d->b[dir2],tree));
2247 PhyML_Fprintf(stderr,"\n. t1: %f t2: %f",tree->times->nd_t[v1->num],tree->times->nd_t[v2->num]);
2248 PhyML_Fprintf(stderr,"\n. rr1: %f rr2: %f",tree->rates->br_r[v1->num],tree->rates->br_r[v2->num]);
2249 PhyML_Fprintf(stderr,"\n. est: %f %f diff: %G",t1,t2,t1-t2);
2250 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
2251 }
2252
2253 tree->times->nd_t[d->num] = t1;
2254 }
2255 }
2256
2257 //////////////////////////////////////////////////////////////
2258 //////////////////////////////////////////////////////////////
2259
TIMES_Copy_Time_Struct(t_time * from,t_time * to,int n_otu)2260 void TIMES_Copy_Time_Struct(t_time *from, t_time *to, int n_otu)
2261 {
2262 int i;
2263
2264 to->birth_rate = from->birth_rate;
2265 to->birth_rate_min = from->birth_rate_min;
2266 to->birth_rate_max = from->birth_rate_max;
2267 to->birth_rate_pivot = from->birth_rate_pivot;
2268
2269 to->death_rate = from->death_rate;
2270 to->death_rate_min = from->death_rate_min;
2271 to->death_rate_max = from->death_rate_max;
2272 to->death_rate_pivot = from->death_rate_pivot;
2273
2274 to->c_lnL_times = from->c_lnL_times;
2275 to->c_lnL_jps = from->c_lnL_jps;
2276
2277 to->nd_t_recorded = from->nd_t_recorded;
2278 to->update_time_norm_const = from->update_time_norm_const;
2279
2280 for(i=0;i<2*n_otu-1;++i) to->nd_t[i] = from->nd_t[i];
2281 for(i=0;i<2*n_otu-1;++i) to->buff_t[i] = from->buff_t[i];
2282 for(i=0;i<2*n_otu-1;++i) to->true_t[i] = from->true_t[i];
2283 for(i=0;i<2*n_otu-1;++i) to->t_mean[i] = from->t_mean[i];
2284 for(i=0;i<2*n_otu-1;++i) to->t_prior[i] = from->t_prior[i];
2285 for(i=0;i<2*n_otu-1;++i) to->t_prior_min[i] = from->t_prior_min[i];
2286 for(i=0;i<2*n_otu-1;++i) to->t_prior_max[i] = from->t_prior_max[i];
2287 for(i=0;i<2*n_otu-1;++i) to->t_floor[i] = from->t_floor[i];
2288 for(i=0;i<2*n_otu-1;++i) to->t_rank[i] = from->t_rank[i];
2289 for(i=0;i<2*n_otu-1;++i) to->t_has_prior[i] = from->t_has_prior[i];
2290 for(i=0;i<2*n_otu-1;++i) to->n_jps[i] = from->n_jps[i];
2291 for(i=0;i<2*n_otu-2;++i) to->t_jps[i] = from->t_jps[i];
2292 for(i=0;i<2*n_otu-1;++i) to->mean_t[i] = from->mean_t[i];
2293 for(i=0;i<2*n_otu-1;++i) to->time_slice_lims[i] = from->time_slice_lims[i];
2294 for(i=0;i<2*n_otu-1;++i) to->n_time_slice_spans[i] = from->n_time_slice_spans[i];
2295 for(i=0;i<2*n_otu-1;++i) to->curr_slice[i] = from->curr_slice[i];
2296 for(i=0;i<2*n_otu-1;++i) to->has_survived[i] = from->has_survived[i];
2297 for(i=0;i<2*n_otu-1;++i) to->calib_prob[i] = from->calib_prob[i];
2298 for(i=0;i<2*n_otu-1;++i) to->t_prior_min_ori[i] = from->t_prior_max_ori[i];
2299 for(i=0;i<n_otu*n_otu;++i) to->times_partial_proba[i] = from->times_partial_proba[i];
2300 for(i=0;i<n_otu*n_otu;++i) to->numb_calib_chosen[i] = from->numb_calib_chosen[i];
2301 for(i=0;i<2*n_otu-1;++i) to->mean_t[i] = from->mean_t[i];
2302 }
2303
2304 //////////////////////////////////////////////////////////////
2305 //////////////////////////////////////////////////////////////
2306
2307 //////////////////////////////////////////////////////////////
2308 //////////////////////////////////////////////////////////////
2309
2310 //////////////////////////////////////////////////////////////
2311 //////////////////////////////////////////////////////////////
2312
2313 //////////////////////////////////////////////////////////////
2314 //////////////////////////////////////////////////////////////
2315
2316 //////////////////////////////////////////////////////////////
2317 //////////////////////////////////////////////////////////////
2318
2319
2320