1 /*
2
3 PhyML: a program that computes maximum likelihood phylogenies from
4 DNA or AA homologous sequences.
5
6 Copyright (C) Stephane Guindon. Oct 2003 onward.
7
8 All parts of the source except where indicated are distributed under
9 the GNU public licence. See http://www.opensource.org for details.
10
11 */
12
13 #include "mixt.h"
14
15 int n_sec1 = 0;
16
17 //////////////////////////////////////////////////////////////
18 //////////////////////////////////////////////////////////////
19
MIXT_Chain_All(t_tree * mixt_tree)20 void MIXT_Chain_All(t_tree *mixt_tree)
21 {
22 t_tree *curr, *next;
23 int i;
24
25 curr = mixt_tree;
26 next = mixt_tree->next;
27
28 do
29 {
30 MIXT_Chain_String(curr->mod->aa_rate_mat_file,next->mod->aa_rate_mat_file);
31 MIXT_Chain_String(curr->mod->modelname,next->mod->modelname);
32 MIXT_Chain_String(curr->mod->custom_mod_string,next->mod->custom_mod_string);
33 MIXT_Chain_Scalar_Dbl(curr->mod->kappa,next->mod->kappa);
34 MIXT_Chain_Scalar_Dbl(curr->mod->lambda,next->mod->lambda);
35 MIXT_Chain_Scalar_Dbl(curr->mod->br_len_mult,next->mod->br_len_mult);
36 MIXT_Chain_Scalar_Dbl(curr->mod->br_len_mult_unscaled,next->mod->br_len_mult_unscaled);
37 MIXT_Chain_Scalar_Dbl(curr->mod->mr,next->mod->mr);
38 MIXT_Chain_Vector_Dbl(curr->mod->Pij_rr,next->mod->Pij_rr);
39 MIXT_Chain_Vector_Dbl(curr->mod->e_frq->user_b_freq,next->mod->e_frq->user_b_freq);
40 for(i=0;i<2*mixt_tree->n_otu-1;++i) MIXT_Chain_Scalar_Dbl(curr->a_edges[i]->l,next->a_edges[i]->l);
41 for(i=0;i<2*mixt_tree->n_otu-1;++i) MIXT_Chain_Scalar_Dbl(curr->a_edges[i]->l_old,next->a_edges[i]->l_old);
42 for(i=0;i<2*mixt_tree->n_otu-1;++i) MIXT_Chain_Scalar_Dbl(curr->a_edges[i]->l_var,next->a_edges[i]->l_var);
43 for(i=0;i<2*mixt_tree->n_otu-1;++i) MIXT_Chain_Scalar_Dbl(curr->a_edges[i]->l_var_old,next->a_edges[i]->l_var_old);
44 MIXT_Chain_Rmat(curr->mod->r_mat,next->mod->r_mat);
45 MIXT_Chain_RAS(curr->mod->ras,next->mod->ras);
46 MIXT_Chain_Efrq(curr->mod->e_frq,next->mod->e_frq);
47 MIXT_Chain_Eigen(curr->mod->eigen,next->mod->eigen);
48
49 curr = next;
50 next = next->next;
51 }
52 while(next);
53
54 MIXT_Chain_Edges(mixt_tree);
55 MIXT_Chain_Nodes(mixt_tree);
56 MIXT_Chain_Sprs(mixt_tree);
57 Make_Rmat_Weight(mixt_tree);
58 Make_Efrq_Weight(mixt_tree);
59
60 }
61
62 //////////////////////////////////////////////////////////////
63 //////////////////////////////////////////////////////////////
64
MIXT_Chain_Edges(t_tree * mixt_tree)65 void MIXT_Chain_Edges(t_tree *mixt_tree)
66 {
67 int i;
68 t_edge *b;
69 t_tree *tree;
70
71 tree = mixt_tree;
72 do
73 {
74 for(i=0;i<2*tree->n_otu-1;++i)
75 {
76 b = tree->a_edges[i];
77
78 if(tree->next) b->next = tree->next->a_edges[i];
79 if(tree->prev) b->prev = tree->prev->a_edges[i];
80 if(tree->next_mixt) b->next_mixt = tree->next_mixt->a_edges[i];
81 if(tree->prev_mixt) b->prev_mixt = tree->prev_mixt->a_edges[i];
82 }
83
84 if(tree->e_root != NULL)
85 {
86 b = tree->e_root;
87
88 if(tree->next) b->next = tree->next->e_root;
89 if(tree->prev) b->prev = tree->prev->e_root;
90 if(tree->next_mixt) b->next_mixt = tree->next_mixt->e_root;
91 if(tree->prev_mixt) b->prev_mixt = tree->prev_mixt->e_root;
92 }
93 tree = tree->next;
94 }
95 while(tree);
96 }
97
98 //////////////////////////////////////////////////////////////
99 //////////////////////////////////////////////////////////////
100
MIXT_Chain_Nodes(t_tree * mixt_tree)101 void MIXT_Chain_Nodes(t_tree *mixt_tree)
102 {
103 int i;
104 t_node *n;
105 t_tree *tree;
106
107 tree = mixt_tree;
108 do
109 {
110 for(i=0;i<2*tree->n_otu-2;++i)
111 {
112 n = tree->a_nodes[i];
113
114 if(tree->next) n->next = tree->next->a_nodes[i];
115 if(tree->prev) n->prev = tree->prev->a_nodes[i];
116 if(tree->next_mixt) n->next_mixt = tree->next_mixt->a_nodes[i];
117 if(tree->prev_mixt) n->prev_mixt = tree->prev_mixt->a_nodes[i];
118 }
119
120 if(tree->n_root != NULL)
121 {
122 n = tree->n_root;
123 if(tree->next) n->next = tree->next->n_root;
124 if(tree->prev) n->prev = tree->prev->n_root;
125 if(tree->next_mixt) n->next_mixt = tree->next_mixt->n_root;
126 if(tree->prev_mixt) n->prev_mixt = tree->prev_mixt->n_root;
127 }
128 tree = tree->next;
129 }
130 while(tree);
131 }
132 //////////////////////////////////////////////////////////////
133 //////////////////////////////////////////////////////////////
134
MIXT_Chain_Sprs(t_tree * mixt_tree)135 void MIXT_Chain_Sprs(t_tree *mixt_tree)
136 {
137 int i;
138 t_tree *tree;
139
140 tree = mixt_tree;
141 do
142 {
143 if(tree->next) tree->best_spr->next = tree->next->best_spr;
144 if(tree->prev) tree->best_spr->prev = tree->prev->best_spr;
145 if(tree->next_mixt) tree->best_spr->next_mixt = tree->next_mixt->best_spr;
146 if(tree->prev_mixt) tree->best_spr->prev_mixt = tree->prev_mixt->best_spr;
147
148 for(i=0;i<2*tree->n_otu-2;++i)
149 {
150 if(tree->next) tree->spr_list_one_edge[i]->next = tree->next->spr_list_one_edge[i];
151 if(tree->prev) tree->spr_list_one_edge[i]->prev = tree->prev->spr_list_one_edge[i];
152 if(tree->next_mixt) tree->spr_list_one_edge[i]->next_mixt = tree->next_mixt->spr_list_one_edge[i];
153 if(tree->prev_mixt) tree->spr_list_one_edge[i]->prev_mixt = tree->prev_mixt->spr_list_one_edge[i];
154
155 if(tree->next) tree->spr_list_all_edge[i]->next = tree->next->spr_list_all_edge[i];
156 if(tree->prev) tree->spr_list_all_edge[i]->prev = tree->prev->spr_list_all_edge[i];
157 if(tree->next_mixt) tree->spr_list_all_edge[i]->next_mixt = tree->next_mixt->spr_list_all_edge[i];
158 if(tree->prev_mixt) tree->spr_list_all_edge[i]->prev_mixt = tree->prev_mixt->spr_list_all_edge[i];
159
160 }
161 tree = tree->next;
162 }
163 while(tree);
164 }
165
166 //////////////////////////////////////////////////////////////
167 //////////////////////////////////////////////////////////////
168
MIXT_Chain_String(t_string * curr,t_string * next)169 void MIXT_Chain_String(t_string *curr, t_string *next)
170 {
171 if(!next)
172 {
173 return;
174 }
175 else
176 {
177 t_string *buff,*last;
178
179 last = NULL;
180
181 /*! Search backward */
182 buff = curr;
183 while(buff)
184 {
185 if(buff == next) break;
186 buff = buff->prev;
187 }
188
189 /*! Search forward */
190 if(!buff)
191 {
192 buff = curr;
193 while(buff)
194 {
195 if(buff == next) break;
196 buff = buff->next;
197 }
198 }
199
200
201 if(!buff)
202 {
203 last = curr;
204 while(last->next) { last = last->next; }
205
206 last->next = next;
207 next->prev = last;
208 }
209 }
210 }
211
212 //////////////////////////////////////////////////////////////
213 //////////////////////////////////////////////////////////////
214
MIXT_Chain_Vector_Dbl(vect_dbl * curr,vect_dbl * next)215 void MIXT_Chain_Vector_Dbl(vect_dbl *curr, vect_dbl *next)
216 {
217 if(!next)
218 {
219 return;
220 }
221 else
222 {
223 vect_dbl *buff,*last;
224
225 last = NULL;
226
227 buff = curr;
228 while(buff)
229 {
230 if(buff == next) break;
231 buff = buff->prev;
232 }
233
234 /*! Search forward */
235 if(!buff)
236 {
237 buff = curr;
238 while(buff)
239 {
240 if(buff == next) break;
241 buff = buff->next;
242 }
243 }
244
245 if(!buff)
246 {
247 last = curr;
248 while(last->next) { last = last->next; }
249
250 last->next = next;
251 next->prev = last;
252 }
253 }
254 }
255
256 //////////////////////////////////////////////////////////////
257 //////////////////////////////////////////////////////////////
258
MIXT_Chain_Scalar_Dbl(scalar_dbl * curr,scalar_dbl * next)259 void MIXT_Chain_Scalar_Dbl(scalar_dbl *curr, scalar_dbl *next)
260 {
261 if(!next)
262 {
263 return;
264 }
265 else
266 {
267 scalar_dbl *buff, *last;
268
269 last = NULL;
270
271 /*! Search backward */
272 buff = curr;
273 while(buff)
274 {
275 if(buff == next) break;
276 buff = buff->prev;
277 }
278
279 /*! Search forward */
280 if(!buff)
281 {
282 buff = curr;
283 while(buff)
284 {
285 if(buff == next) break;
286 buff = buff->next;
287 }
288 }
289
290 /*! Not chained yet. Add next at the end of chained list */
291 if(!buff)
292 {
293 last = curr;
294 while(last->next) { last = last->next; }
295
296 last->next = next;
297 next->prev = last;
298 }
299 }
300 }
301
302 //////////////////////////////////////////////////////////////
303 //////////////////////////////////////////////////////////////
304
MIXT_Chain_Rmat(t_rmat * curr,t_rmat * next)305 void MIXT_Chain_Rmat(t_rmat *curr, t_rmat *next)
306 {
307 if(!next)
308 {
309 return;
310 }
311 else
312 {
313 t_rmat *buff, *last;
314
315 last = NULL;
316
317 buff = curr;
318 while(buff)
319 {
320 if(buff == next) break;
321 buff = buff->prev;
322 }
323
324 /*! Search forward */
325 if(!buff)
326 {
327 buff = curr;
328 while(buff)
329 {
330 if(buff == next) break;
331 buff = buff->next;
332 }
333 }
334
335 if(!buff)
336 {
337 last = curr;
338 while(last->next) { last = last->next; }
339
340 last->next = next;
341 next->prev = last;
342 }
343 }
344 }
345
346 //////////////////////////////////////////////////////////////
347 //////////////////////////////////////////////////////////////
348
MIXT_Chain_Efrq(t_efrq * curr,t_efrq * next)349 void MIXT_Chain_Efrq(t_efrq *curr, t_efrq *next)
350 {
351 if(!next)
352 {
353 return;
354 }
355 else
356 {
357 t_efrq *buff,*last;
358
359 last = NULL;
360
361 buff = curr;
362 while(buff)
363 {
364 if(buff == next) break;
365 buff = buff->prev;
366 }
367
368 /*! Search forward */
369 if(!buff)
370 {
371 buff = curr;
372 while(buff)
373 {
374 if(buff == next) break;
375 buff = buff->next;
376 }
377 }
378
379 if(!buff)
380 {
381 last = curr;
382 while(last->next) { last = last->next; }
383
384 last->next = next;
385 next->prev = last;
386 }
387 }
388 }
389
390 //////////////////////////////////////////////////////////////
391 //////////////////////////////////////////////////////////////
392
MIXT_Chain_Eigen(eigen * curr,eigen * next)393 void MIXT_Chain_Eigen(eigen *curr, eigen *next)
394 {
395 if(!next)
396 {
397 return;
398 }
399 else
400 {
401 eigen *buff,*last;
402
403 last = NULL;
404
405 buff = curr;
406 while(buff)
407 {
408 if(buff == next) break;
409 buff = buff->prev;
410 }
411
412 /*! Search forward */
413 if(!buff)
414 {
415 buff = curr;
416 while(buff)
417 {
418 if(buff == next) break;
419 buff = buff->next;
420 }
421 }
422
423 if(!buff)
424 {
425 last = curr;
426 while(last->next) { last = last->next; }
427
428 last->next = next;
429 next->prev = last;
430 }
431 }
432 }
433
434 //////////////////////////////////////////////////////////////
435 //////////////////////////////////////////////////////////////
436
MIXT_Chain_RAS(t_ras * curr,t_ras * next)437 void MIXT_Chain_RAS(t_ras *curr, t_ras *next)
438 {
439 if(!next) return;
440 else
441 {
442 t_ras *buff,*last;
443
444 last = NULL;
445
446 buff = curr;
447 while(buff)
448 {
449 if(buff == next) break;
450 buff = buff->prev;
451 }
452
453 /*! Search forward */
454 if(!buff)
455 {
456 buff = curr;
457 while(buff)
458 {
459 if(buff == next) break;
460 buff = buff->next;
461 }
462 }
463
464 if(!buff)
465 {
466 last = curr;
467 while(last->next) { last = last->next; }
468
469 last->next = next;
470 next->prev = last;
471 }
472 }
473 }
474
475 //////////////////////////////////////////////////////////////
476 //////////////////////////////////////////////////////////////
477
MIXT_Chain_Rates(t_rate * curr,t_rate * next)478 void MIXT_Chain_Rates(t_rate *curr, t_rate *next)
479 {
480 if(!next) return;
481 else
482 {
483 t_rate *buff,*last;
484
485 last = NULL;
486
487 buff = curr;
488 while(buff)
489 {
490 if(buff == next) break;
491 buff = buff->prev;
492 }
493
494 /*! Search forward */
495 if(!buff)
496 {
497 buff = curr;
498 while(buff)
499 {
500 if(buff == next) break;
501 buff = buff->next;
502 }
503 }
504
505 if(!buff)
506 {
507 last = curr;
508 while(last->next) { last = last->next; }
509
510 last->next = next;
511 next->prev = last;
512 }
513 }
514 }
515
516 //////////////////////////////////////////////////////////////
517 //////////////////////////////////////////////////////////////
518
MIXT_Chain_Cal(t_tree * mixt_tree)519 void MIXT_Chain_Cal(t_tree *mixt_tree)
520 {
521 int i;
522
523 for(i=0;i<mixt_tree->times->n_cal-1;i++)
524 {
525 mixt_tree->times->a_cal[i]->next = mixt_tree->times->a_cal[i+1];
526 mixt_tree->times->a_cal[i+1]->prev = mixt_tree->times->a_cal[i];
527 }
528 }
529
530 //////////////////////////////////////////////////////////////
531 //////////////////////////////////////////////////////////////
532
MIXT_Turn_Branches_OnOff_In_All_Elem(int onoff,t_tree * mixt_tree)533 void MIXT_Turn_Branches_OnOff_In_All_Elem(int onoff, t_tree *mixt_tree)
534 {
535 t_tree *tree;
536
537 /*! Turn all branches to ON state */
538 tree = mixt_tree;
539 do
540 {
541 MIXT_Turn_Branches_OnOff_In_One_Elem(ON,tree);
542 tree = tree->next_mixt;
543 }
544 while(tree);
545 }
546
547 //////////////////////////////////////////////////////////////
548 //////////////////////////////////////////////////////////////
549
MIXT_Turn_Branches_OnOff_In_One_Elem(int onoff,t_tree * mixt_tree)550 void MIXT_Turn_Branches_OnOff_In_One_Elem(int onoff, t_tree *mixt_tree)
551 {
552 int i;
553 t_tree *tree;
554
555 if(mixt_tree->is_mixt_tree == NO)
556 {
557 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
558 Exit("\n");
559 }
560
561 tree = mixt_tree;
562
563 do
564 {
565 for(i=0;i<2*tree->n_otu-1;++i) tree->a_edges[i]->l->onoff = onoff;
566 tree = tree->next;
567 }
568 while(tree && tree->is_mixt_tree == NO);
569 }
570
571 //////////////////////////////////////////////////////////////
572 //////////////////////////////////////////////////////////////
573
MIXT_Post_Order_Lk(t_node * mixt_a,t_node * mixt_d,t_tree * mixt_tree)574 void MIXT_Post_Order_Lk(t_node *mixt_a, t_node *mixt_d, t_tree *mixt_tree)
575 {
576 t_tree *tree;
577 t_node *a,*d;
578
579
580 tree = mixt_tree;
581 a = mixt_a;
582 d = mixt_d;
583
584 assert(a);
585 assert(d);
586 assert(tree);
587
588 do
589 {
590 if(tree->is_mixt_tree)
591 {
592 tree = tree->next;
593 a = a->next;
594 d = d->next;
595 }
596
597 assert(a);
598 assert(d);
599 assert(tree);
600
601 if(tree->mod->ras->invar == NO) Post_Order_Lk(a,d,tree);
602
603 tree = tree->next;
604 a = a->next;
605 d = d->next;
606 }
607 while(tree);
608
609 }
610
611 //////////////////////////////////////////////////////////////
612 //////////////////////////////////////////////////////////////
613
MIXT_Pre_Order_Lk(t_node * mixt_a,t_node * mixt_d,t_tree * mixt_tree)614 void MIXT_Pre_Order_Lk(t_node *mixt_a, t_node *mixt_d, t_tree *mixt_tree)
615 {
616 t_tree *tree;
617 t_node *a,*d;
618
619 tree = mixt_tree;
620 a = mixt_a;
621 d = mixt_d;
622
623 assert(a);
624 assert(d);
625 assert(tree);
626
627 do
628 {
629 if(tree->is_mixt_tree)
630 {
631 tree = tree->next;
632 a = a->next;
633 d = d->next;
634 }
635
636
637 assert(a);
638 assert(d);
639 assert(tree);
640
641 if(tree->mod->ras->invar == NO) Pre_Order_Lk(a,d,tree);
642
643 tree = tree->next;
644 a = a->next;
645 d = d->next;
646 }
647 while(tree);
648
649 }
650
651 //////////////////////////////////////////////////////////////
652 //////////////////////////////////////////////////////////////
653
MIXT_Lk(t_edge * mixt_b,t_tree * mixt_tree)654 phydbl MIXT_Lk(t_edge *mixt_b, t_tree *mixt_tree)
655 {
656 t_tree *tree,*cpy_mixt_tree;
657 t_edge *b,*cpy_mixt_b;
658 phydbl sum_lnL;
659 unsigned int site, br, ns;
660 phydbl *sum_scale_left_cat,*sum_scale_rght_cat;
661 phydbl sum;
662 phydbl site_lk_cat,site_lk,log_site_lk,inv_site_lk;
663 int num_prec_issue;
664 int ambiguity_check,state;
665 int l;
666 phydbl r_mat_weight_sum, e_frq_weight_sum, sum_probas;
667 phydbl len;
668 phydbl *expl,*dot_prod;
669
670 tree = NULL;
671 b = NULL;
672 expl = NULL;
673 cpy_mixt_tree = mixt_tree;
674 cpy_mixt_b = mixt_b;
675 len = -1.;
676 ns = 0;
677
678 /* MIXT_Update_Br_Len_Multipliers(mixt_tree->mod); */
679
680 #if (defined PHYTIME || defined INVITEE || defined PHYREX || defined DATE)
681 if((mixt_tree->rates) && (mixt_tree->rates->bl_from_rt)) RATES_Update_Cur_Bl(mixt_tree);
682 #endif
683
684
685 /* Update RAS structure (mixt_tree level) */
686 tree = mixt_tree;
687 do
688 {
689 Update_RAS(tree->mod);
690 tree = tree->next_mixt;
691 }
692 while(tree);
693
694 /* Update other model structure (tree level) */
695 tree = mixt_tree->next;
696 do
697 {
698 if(tree->is_mixt_tree == YES) tree = tree->next;
699
700 if(!cpy_mixt_b)
701 {
702 if(!Update_Boundaries(tree->mod) ||
703 !Update_Efrq(tree->mod) ||
704 !Update_Eigen(tree->mod))
705 {
706 PhyML_Fprintf(stderr,"\n. move: %s",mixt_tree->mcmc->move_name[mixt_tree->mcmc->move_idx]);
707 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
708 }
709 }
710
711 tree = tree->next;
712 }
713 while(tree);
714
715
716 if(!cpy_mixt_b)
717 {
718 for(br=0;br<2*mixt_tree->n_otu-3;++br) MIXT_Update_PMat_At_Given_Edge(mixt_tree->a_edges[br],mixt_tree);
719
720 if(mixt_tree->n_root && mixt_tree->ignore_root == NO)
721 {
722 MIXT_Update_PMat_At_Given_Edge(mixt_tree->n_root->b[1],mixt_tree);
723 MIXT_Update_PMat_At_Given_Edge(mixt_tree->n_root->b[2],mixt_tree);
724 }
725 }
726 else
727 {
728 MIXT_Update_PMat_At_Given_Edge(mixt_b,mixt_tree);
729 }
730
731 if(!cpy_mixt_b)
732 {
733 if(mixt_tree->n_root != NULL)
734 {
735 if(mixt_tree->ignore_root == NO)
736 {
737 MIXT_Post_Order_Lk(mixt_tree->n_root,mixt_tree->n_root->v[1],mixt_tree);
738 MIXT_Post_Order_Lk(mixt_tree->n_root,mixt_tree->n_root->v[2],mixt_tree);
739
740 MIXT_Update_Partial_Lk(mixt_tree,mixt_tree->n_root->b[1],mixt_tree->n_root);
741 MIXT_Update_Partial_Lk(mixt_tree,mixt_tree->n_root->b[2],mixt_tree->n_root);
742
743 if(mixt_tree->both_sides == YES)
744 {
745 MIXT_Pre_Order_Lk(mixt_tree->n_root,mixt_tree->n_root->v[1],mixt_tree);
746 MIXT_Pre_Order_Lk(mixt_tree->n_root,mixt_tree->n_root->v[2],mixt_tree);
747 }
748 }
749 else
750 {
751 MIXT_Post_Order_Lk(mixt_tree->e_root->rght,mixt_tree->e_root->left,mixt_tree);
752 MIXT_Post_Order_Lk(mixt_tree->e_root->left,mixt_tree->e_root->rght,mixt_tree);
753
754 if(mixt_tree->both_sides == YES)
755 {
756 MIXT_Pre_Order_Lk(mixt_tree->e_root->rght,mixt_tree->e_root->left,mixt_tree);
757 MIXT_Pre_Order_Lk(mixt_tree->e_root->left,mixt_tree->e_root->rght,mixt_tree);
758 }
759 }
760 }
761 else
762 {
763 MIXT_Post_Order_Lk(mixt_tree->a_nodes[0],mixt_tree->a_nodes[0]->v[0],mixt_tree);
764 if(mixt_tree->both_sides == YES)
765 {
766 MIXT_Pre_Order_Lk(mixt_tree->a_nodes[0],mixt_tree->a_nodes[0]->v[0],mixt_tree);
767 }
768 }
769 }
770
771 do /*! Consider each element of the data partition */
772 {
773 mixt_tree->c_lnL = 0.0;
774 tree = mixt_tree->next;
775 do
776 {
777 tree->c_lnL = 0.0;
778 tree = tree->next;
779 }
780 while(tree && tree->is_mixt_tree == NO);
781
782 Set_Br_Len_Var(mixt_b,mixt_tree);
783
784 if(!cpy_mixt_b)
785 {
786 if(mixt_tree->n_root)
787 {
788 if(mixt_tree->ignore_root == NO)
789 mixt_b = (mixt_tree->n_root->v[1]->tax == NO)?(mixt_tree->n_root->b[2]):(mixt_tree->n_root->b[1]);
790 else
791 mixt_b = mixt_tree->e_root;
792 }
793 else
794 {
795 mixt_b = mixt_tree->a_nodes[0]->b[0];
796 }
797 }
798
799 sum_scale_left_cat = (phydbl *)mCalloc(MAX(mixt_tree->mod->ras->n_catg,mixt_tree->mod->n_mixt_classes),sizeof(phydbl));
800 sum_scale_rght_cat = (phydbl *)mCalloc(MAX(mixt_tree->mod->ras->n_catg,mixt_tree->mod->n_mixt_classes),sizeof(phydbl));
801
802 r_mat_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(mixt_tree->next->mod->r_mat_weight);
803 e_frq_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(mixt_tree->next->mod->e_frq_weight);
804 sum_probas = MIXT_Get_Sum_Of_Probas_Across_Mixtures(r_mat_weight_sum,e_frq_weight_sum,mixt_tree);
805
806
807 b = mixt_b->next;
808 tree = mixt_tree->next;
809 do
810 {
811 while(tree->mod->ras->invar == YES)
812 {
813 tree = tree->next;
814 b = b->next;
815
816 if(tree == NULL || tree->is_mixt_tree == YES)
817 {
818 PhyML_Fprintf(stderr,"\n. %p",(void *)tree);
819 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
820 Exit("\n");
821 }
822 }
823
824 if(tree->update_eigen_lr == YES) Update_Eigen_Lr(b,tree);
825
826 if(tree->use_eigen_lr == YES)
827 {
828 len = b->l->v;
829 len *= tree->mod->br_len_mult->v;
830 len *= tree->mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number];
831 if(len < tree->mod->l_min) len = tree->mod->l_min;
832 else if(len > tree->mod->l_max) len = tree->mod->l_max;
833 expl = tree->expl;
834 for(l=0;l<tree->mod->ns;l++) expl[l] = (phydbl)POW(tree->mod->eigen->e_val[l],len);
835 }
836
837 tree = tree->next;
838 b = b->next;
839
840 }
841 while(tree && tree->is_mixt_tree == NO);
842
843 tree = mixt_tree;
844 do
845 {
846 tree->numerical_warning = NO;
847 tree = tree->next;
848 }
849 while(tree);
850
851 mixt_tree->c_lnL = .0;
852
853 for(site=0;site<mixt_tree->n_pattern;++site)
854 {
855 b = mixt_b->next;
856 tree = mixt_tree->next;
857
858 /*! Skip calculations if model has zero rate */
859 while(tree->mod->ras->invar == YES)
860 {
861 tree = tree->next;
862 b = b->next;
863
864 if(tree == NULL || tree->is_mixt_tree == YES)
865 {
866 PhyML_Fprintf(stderr,"\n. %p",(void *)tree);
867 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
868 Exit("\n");
869 }
870 }
871
872 ambiguity_check = -1;
873 state = -1;
874
875 if((b->rght->tax) &&
876 (!mixt_tree->mod->s_opt->greedy) &&
877 (mixt_tree->data->wght[site] > SMALL))
878 {
879 ambiguity_check = b->rght->c_seq->is_ambigu[site];
880 if(!ambiguity_check) state = b->rght->c_seq->d_state[site];
881 }
882
883 /*! For all classes in the mixture */
884 do
885 {
886 if(tree->is_mixt_tree == YES)
887 {
888 tree = tree->next;
889 b = b->next;
890 }
891
892 ns = tree->mod->ns;
893 tree->curr_site = site;
894 tree->apply_lk_scaling = NO;
895 site_lk_cat = 0.0;
896 expl = tree->expl;
897 dot_prod = tree->dot_prod;
898
899 if(!(tree->mod->ras->invar == YES && mixt_tree->is_mixt_tree == YES) && (tree->data->wght[tree->curr_site] > SMALL))
900 {
901 if((b->rght->tax) && (tree->mod->s_opt->greedy == NO))
902 {
903 ambiguity_check = b->rght->c_seq->is_ambigu[tree->curr_site];
904 if(ambiguity_check == NO) state = b->rght->c_seq->d_state[tree->curr_site];
905 }
906
907 if(tree->use_eigen_lr == YES)
908 {
909 site_lk_cat = Lk_Core_Eigen_Lr(expl,dot_prod + site*ns,b,tree);
910 }
911 else
912 {
913 if(b->rght->tax == YES)
914 site_lk_cat = Lk_Core(state,ambiguity_check,
915 b->p_lk_left + site*ns,
916 b->p_lk_tip_r + site*ns,
917 b->Pij_rr,b->tPij_rr,b,tree);
918 else
919 site_lk_cat = Lk_Core(state,ambiguity_check,
920 b->p_lk_left + site*ns,
921 b->p_lk_rght + site*ns,
922 b->Pij_rr,b->tPij_rr,b,tree);
923 }
924 }
925
926
927 tree->apply_lk_scaling = YES;
928
929 sum_scale_left_cat[tree->mod->ras->parent_class_number] =
930 (b->sum_scale_left)?
931 (b->sum_scale_left[site]):
932 (0.0);
933
934 sum_scale_rght_cat[tree->mod->ras->parent_class_number] =
935 (b->sum_scale_rght)?
936 (b->sum_scale_rght[site]):
937 (0.0);
938
939 sum =
940 sum_scale_left_cat[tree->mod->ras->parent_class_number] +
941 sum_scale_rght_cat[tree->mod->ras->parent_class_number];
942
943 if(sum > 1024.)
944 {
945 /* PhyML_Fprintf(stderr,"\n. Numerical precision issue detected (sum = %g)!!!",sum); */
946 sum = 1023.;
947 tree->mixt_tree->numerical_warning = YES;
948 tree->numerical_warning = YES;
949 }
950
951 tree->unscaled_site_lk_cat[site] = site_lk_cat;
952
953 site_lk_cat /= pow(2,sum);
954
955 tree->site_lk_cat[0] = site_lk_cat;
956
957 tree = tree->next;
958 b = b->next;
959 }
960 while(tree && tree->is_mixt_tree == NO);
961 // done with all trees in the mixture for this partition element.
962 // All likelihood values are in site_lk_cat[0]
963
964
965 tree = mixt_tree->next;
966 b = mixt_b->next;
967 site_lk = .0;
968
969 do
970 {
971 if(tree->mod->ras->invar == YES)
972 {
973 tree = tree->next;
974 b = b->next;
975 if(!(tree && tree->is_mixt_tree == NO)) break;
976 }
977
978 site_lk +=
979 tree->site_lk_cat[0] *
980 mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number] *
981 tree->mod->r_mat_weight->v / r_mat_weight_sum *
982 tree->mod->e_frq_weight->v / e_frq_weight_sum /
983 sum_probas;
984
985 tree = tree->next;
986 b = b->next;
987 }
988 while(tree && tree->is_mixt_tree == NO);
989
990
991 /* Scaling for invariants */
992 if(mixt_tree->mod->ras->invar == YES)
993 {
994 num_prec_issue = NO;
995
996 tree = mixt_tree->next;
997 while(tree->mod->ras->invar == NO)
998 {
999 tree = tree->next;
1000 if(tree == NULL || tree->is_mixt_tree == YES)
1001 {
1002 PhyML_Fprintf(stderr,"\n. tree: %p",tree);
1003 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d",__FILE__,__LINE__);
1004 Exit("\n");
1005 }
1006 }
1007
1008 /*! 'tree' will give the correct state frequencies (as opposed to mixt_tree) */
1009 inv_site_lk = Invariant_Lk(0,site,&num_prec_issue,tree);
1010
1011 if(num_prec_issue == YES) // inv_site_lk >> site_lk
1012 {
1013 site_lk = inv_site_lk * mixt_tree->mod->ras->pinvar->v;
1014 }
1015 else
1016 {
1017 site_lk = site_lk * (1. - mixt_tree->mod->ras->pinvar->v) + inv_site_lk * mixt_tree->mod->ras->pinvar->v;
1018 }
1019 }
1020
1021 if(site_lk < SMALL)
1022 {
1023 /* PhyML_Fprintf(stderr,"\n. site = %d",site); */
1024 /* PhyML_Fprintf(stderr,"\n. invar = %d",mixt_tree->data->invar[site]); */
1025 /* PhyML_Fprintf(stderr,"\n. mixt = %d",mixt_tree->is_mixt_tree); */
1026 /* PhyML_Fprintf(stderr,"\n. lk = %G log(lk) = %f < %G",site_lk,log_site_lk,-BIG); */
1027 /* for(class=0;class<mixt_tree->mod->ras->n_catg;class++) PhyML_Fprintf(stderr,"\n. rr=%f p=%f",mixt_tree->mod->ras->gamma_rr->v[class],mixt_tree->mod->ras->gamma_r_proba->v[class]); */
1028 /* PhyML_Fprintf(stderr,"\n. pinv = %G",mixt_tree->mod->ras->pinvar->v); */
1029 /* PhyML_Fprintf(stderr,"\n. bl mult = %G",mixt_tree->mod->br_len_mult->v); */
1030 /* PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d.\n",__FILE__,__LINE__); */
1031 /* Exit("\n"); */
1032
1033 site_lk = SMALL;
1034 mixt_tree->numerical_warning = YES;
1035 }
1036
1037 log_site_lk = log(site_lk);
1038
1039
1040
1041 // ... or using the log-likelihood
1042 if(isinf(site_lk) || isnan(site_lk))
1043 {
1044 mixt_tree->cur_site_lk[site] = exp(log_site_lk);
1045 }
1046 else
1047 {
1048 mixt_tree->cur_site_lk[site] = site_lk;
1049 }
1050
1051
1052 /* Multiply log likelihood by the number of times this site pattern is found in the data */
1053 mixt_tree->c_lnL_sorted[site] = mixt_tree->data->wght[site]*log_site_lk;
1054
1055
1056 mixt_tree->c_lnL += mixt_tree->data->wght[site]*log_site_lk;
1057
1058 }
1059
1060 Free(sum_scale_left_cat);
1061 Free(sum_scale_rght_cat);
1062
1063 mixt_tree = mixt_tree->next_mixt;
1064 mixt_b = mixt_b->next_mixt;
1065 }
1066 while(mixt_tree);
1067
1068 mixt_tree = cpy_mixt_tree;
1069 mixt_b = cpy_mixt_b;
1070
1071 sum_lnL = .0;
1072 do
1073 {
1074 sum_lnL += mixt_tree->c_lnL;
1075 mixt_tree = mixt_tree->next_mixt;
1076 }
1077 while(mixt_tree);
1078
1079 mixt_tree = cpy_mixt_tree;
1080 do
1081 {
1082 mixt_tree->c_lnL = sum_lnL;
1083 mixt_tree = mixt_tree->next_mixt;
1084 }
1085 while(mixt_tree);
1086
1087 mixt_tree = cpy_mixt_tree;
1088
1089 return mixt_tree->c_lnL;
1090 }
1091
1092 //////////////////////////////////////////////////////////////
1093 //////////////////////////////////////////////////////////////
1094
MIXT_Update_Partial_Lk(t_tree * mixt_tree,t_edge * mixt_b,t_node * mixt_d)1095 void MIXT_Update_Partial_Lk(t_tree *mixt_tree, t_edge *mixt_b, t_node *mixt_d)
1096 {
1097 t_tree *tree;
1098 t_edge *b;
1099 t_node *d;
1100
1101 tree = mixt_tree;
1102 b = mixt_b;
1103 d = mixt_d;
1104
1105 do
1106 {
1107 if(tree->is_mixt_tree)
1108 {
1109 tree = tree->next;
1110 b = b->next;
1111 d = d->next;
1112 }
1113
1114 if(tree->mod->ras->invar == NO) Update_Partial_Lk(tree,b,d);
1115
1116 tree = tree->next;
1117 b = b->next;
1118 d = d->next;
1119
1120 }
1121 while(tree);
1122
1123 }
1124
1125 //////////////////////////////////////////////////////////////
1126 //////////////////////////////////////////////////////////////
1127
MIXT_Update_Partial_Pars(t_tree * mixt_tree,t_edge * mixt_b,t_node * mixt_d)1128 void MIXT_Update_Partial_Pars(t_tree *mixt_tree, t_edge *mixt_b, t_node *mixt_d)
1129 {
1130 t_tree *tree;
1131 t_edge *b;
1132 t_node *d;
1133
1134 tree = mixt_tree;
1135 b = mixt_b;
1136 d = mixt_d;
1137
1138 do
1139 {
1140 if(tree->is_mixt_tree)
1141 {
1142 tree = tree->next;
1143 b = b->next;
1144 d = d->next;
1145 }
1146
1147 Update_Partial_Pars(tree,b,d);
1148
1149 tree = tree->next;
1150 b = b->next;
1151 d = d->next;
1152
1153 }
1154 while(tree);
1155
1156 }
1157
1158 //////////////////////////////////////////////////////////////
1159 //////////////////////////////////////////////////////////////
1160
MIXT_Update_PMat_At_Given_Edge(t_edge * mixt_b,t_tree * mixt_tree)1161 void MIXT_Update_PMat_At_Given_Edge(t_edge *mixt_b, t_tree *mixt_tree)
1162 {
1163 t_tree *tree;
1164 t_edge *b;
1165
1166 tree = mixt_tree;
1167 b = mixt_b;
1168
1169 do
1170 {
1171 if(tree->is_mixt_tree)
1172 {
1173 tree = tree->next;
1174 b = b->next;
1175 }
1176
1177 if(tree->mod->ras->invar == NO) Update_PMat_At_Given_Edge(b,tree);
1178
1179 tree = tree->next;
1180 b = b->next;
1181 }
1182 while(tree);
1183 }
1184
1185 //////////////////////////////////////////////////////////////
1186 //////////////////////////////////////////////////////////////
1187
MIXT_Get_Number_Of_Classes_In_All_Mixtures(t_tree * mixt_tree)1188 int *MIXT_Get_Number_Of_Classes_In_All_Mixtures(t_tree *mixt_tree)
1189 {
1190 int *n_catg;
1191 t_tree *tree;
1192 int class;
1193
1194 if(mixt_tree->is_mixt_tree == YES)
1195 {
1196 n_catg = NULL;
1197 tree = mixt_tree;
1198 class = 0;
1199 do
1200 {
1201 if(!class) n_catg = (int *)mCalloc(1,sizeof(int));
1202 else n_catg = (int *)realloc(n_catg,(class+1)*sizeof(int));
1203
1204 tree = tree->next;
1205
1206 n_catg[class]=0;
1207 do
1208 {
1209 n_catg[class]++;
1210 tree = tree->next;
1211 }
1212 while(tree && tree->is_mixt_tree == NO);
1213
1214 class++;
1215 }
1216 while(tree);
1217 }
1218 else
1219 {
1220 n_catg = (int *)mCalloc(1,sizeof(int));
1221 n_catg[0] = mixt_tree->mod->ras->n_catg;
1222 if(mixt_tree->mod->ras->invar == YES) n_catg[0]++;
1223 }
1224 return(n_catg);
1225 }
1226
1227 //////////////////////////////////////////////////////////////
1228 //////////////////////////////////////////////////////////////
1229
MIXT_Record_All_Mixtures(t_tree * mixt_tree)1230 t_tree **MIXT_Record_All_Mixtures(t_tree *mixt_tree)
1231 {
1232 t_tree **tree_list;
1233 int n_trees;
1234 t_tree *tree;
1235
1236 tree_list = NULL;
1237 n_trees = 0;
1238 tree = mixt_tree;
1239 do
1240 {
1241 if(!tree_list) tree_list = (t_tree **)mCalloc(1,sizeof(t_tree *));
1242 else tree_list = (t_tree **)realloc(tree_list,(n_trees+1)*sizeof(t_tree *));
1243
1244 tree_list[n_trees] = tree;
1245 n_trees++;
1246 tree = tree->next;
1247 }
1248 while(tree);
1249
1250 tree_list = (t_tree **)realloc(tree_list,(n_trees+1)*sizeof(t_tree *));
1251 tree_list[n_trees] = NULL;
1252
1253 return(tree_list);
1254 }
1255
1256 //////////////////////////////////////////////////////////////
1257 //////////////////////////////////////////////////////////////
1258
MIXT_Break_All_Mixtures(int * c_max,t_tree * mixt_tree)1259 void MIXT_Break_All_Mixtures(int *c_max, t_tree *mixt_tree)
1260 {
1261 t_tree *tree;
1262 int c,i,n;
1263
1264 if(mixt_tree->is_mixt_tree == NO) return;
1265
1266 c = 0;
1267 n = -1;
1268 tree = mixt_tree;
1269 do
1270 {
1271 if(tree->is_mixt_tree == YES)
1272 {
1273 c = 0;
1274 n++;
1275 tree = tree->next;
1276 }
1277
1278 if(c == (c_max[n]-1) &&
1279 tree->next != NULL &&
1280 tree->next->is_mixt_tree == NO)
1281 {
1282 if(tree->mixt_tree->next_mixt == NULL)
1283 {
1284 tree->next = NULL;
1285 for(i=0;i<2*tree->n_otu-1;++i) tree->a_edges[i]->next = NULL;
1286 for(i=0;i<2*tree->n_otu-1;++i) tree->a_nodes[i]->next = NULL;
1287 for(i=0;i<2*tree->n_otu-2;++i) tree->spr_list_one_edge[i]->next = NULL;
1288 }
1289 else
1290 {
1291 tree->next = tree->mixt_tree->next_mixt;
1292 for(i=0;i<2*tree->n_otu-1;++i) tree->a_edges[i]->next = tree->mixt_tree->next_mixt->a_edges[i];
1293 for(i=0;i<2*tree->n_otu-1;++i) tree->a_nodes[i]->next = tree->mixt_tree->next_mixt->a_nodes[i];
1294 for(i=0;i<2*tree->n_otu-2;++i) tree->spr_list_one_edge[i]->next = tree->mixt_tree->next_mixt->spr_list_one_edge[i];
1295 }
1296 }
1297
1298 tree = tree->next;
1299 c++;
1300 }
1301 while(tree);
1302 }
1303
1304 //////////////////////////////////////////////////////////////
1305 //////////////////////////////////////////////////////////////
1306
MIXT_Reconnect_All_Mixtures(t_tree ** tree_list,t_tree * mixt_tree)1307 void MIXT_Reconnect_All_Mixtures(t_tree **tree_list, t_tree *mixt_tree)
1308 {
1309 t_tree *tree;
1310 int n_trees;
1311
1312 if(mixt_tree->is_mixt_tree == NO) return;
1313
1314 tree = mixt_tree;
1315 n_trees = 0;
1316 do
1317 {
1318 tree = tree_list[n_trees];
1319 if(tree->is_mixt_tree == NO) tree->next = tree_list[n_trees+1];
1320 n_trees++;
1321 tree = tree->next;
1322 }
1323 while(tree);
1324
1325 MIXT_Chain_All(mixt_tree);
1326 }
1327
1328 //////////////////////////////////////////////////////////////
1329 //////////////////////////////////////////////////////////////
1330
MIXT_Record_Has_Invariants(t_tree * mixt_tree)1331 int *MIXT_Record_Has_Invariants(t_tree *mixt_tree)
1332 {
1333 int *has_invariants;
1334 t_tree *tree;
1335 int n_trees;
1336
1337 has_invariants = NULL;
1338 tree = mixt_tree;
1339 n_trees = 0;
1340 do
1341 {
1342 if(!n_trees) has_invariants = (int *)mCalloc(1,sizeof(int));
1343 else has_invariants = (int *)realloc(has_invariants,(n_trees+1)*sizeof(int));
1344 has_invariants[n_trees] = (tree->mod->ras->invar == YES)?1:0;
1345 n_trees++;
1346 tree = tree->next;
1347 }
1348 while(tree);
1349
1350 return(has_invariants);
1351 }
1352
1353 //////////////////////////////////////////////////////////////
1354 //////////////////////////////////////////////////////////////
1355
MIXT_Reset_Has_Invariants(int * has_invariants,t_tree * mixt_tree)1356 void MIXT_Reset_Has_Invariants(int *has_invariants, t_tree *mixt_tree)
1357 {
1358 t_tree *tree;
1359 int n_trees;
1360
1361 tree = mixt_tree;
1362 n_trees = 0;
1363 do
1364 {
1365 tree->mod->ras->invar = has_invariants[n_trees];
1366 n_trees++;
1367 tree = tree->next;
1368 }
1369 while(tree);
1370
1371 }
1372
1373 //////////////////////////////////////////////////////////////
1374 //////////////////////////////////////////////////////////////
1375
MIXT_Check_Invar_Struct_In_Each_Partition_Elem(t_tree * mixt_tree)1376 void MIXT_Check_Invar_Struct_In_Each_Partition_Elem(t_tree *mixt_tree)
1377 {
1378 if(mixt_tree->is_mixt_tree == NO) return;
1379 else
1380 {
1381 t_tree *tree;
1382 int n_inv;
1383
1384 n_inv = 0;
1385 tree = mixt_tree;
1386 do
1387 {
1388 if(tree->is_mixt_tree)
1389 {
1390 tree = tree->next;
1391 n_inv = 0;
1392 }
1393
1394 if(tree->mod->ras->invar == YES) n_inv++;
1395
1396 if(n_inv > 1)
1397 {
1398 PhyML_Fprintf(stderr,"\n. Found %d classes of the mixture for file '%s' set to",n_inv,tree->mixt_tree->io->in_align_file);
1399 PhyML_Fprintf(stderr,"\n. invariable. Only one such class per mixture is allowed.");
1400 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
1401 Warn_And_Exit("\n");
1402 }
1403
1404 if(tree->mixt_tree->mod->ras->invar == NO &&
1405 tree->mod->ras->invar == YES)
1406 {
1407 PhyML_Fprintf(stderr,"\n. Unexpected settings for 'siterates' in a partition element (file '%s')",tree->mixt_tree->io->in_align_file);
1408 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
1409 Warn_And_Exit("\n");
1410 }
1411
1412 tree = tree->next;
1413 }
1414 while(tree);
1415 }
1416 }
1417
1418 //////////////////////////////////////////////////////////////
1419 //////////////////////////////////////////////////////////////
1420
MIXT_Check_RAS_Struct_In_Each_Partition_Elem(t_tree * mixt_tree)1421 void MIXT_Check_RAS_Struct_In_Each_Partition_Elem(t_tree *mixt_tree)
1422 {
1423 if(mixt_tree->is_mixt_tree == NO) return;
1424 else
1425 {
1426 t_tree *tree;
1427 int n_classes;
1428
1429 n_classes = 0;
1430 tree = mixt_tree;
1431 do
1432 {
1433 if(tree->is_mixt_tree)
1434 {
1435 if(tree->mod->ras->invar == YES)
1436 {
1437 if(tree->next->mod->ras->invar == NO)
1438 {
1439 PhyML_Fprintf(stderr,"\n. The invariant site class has to be the first element in");
1440 PhyML_Fprintf(stderr,"\n. each <mixtureelem> component. Please amend you XML");
1441 PhyML_Fprintf(stderr,"\n. file accordingly.\n");
1442 Exit("\n.");
1443 }
1444 }
1445 tree = tree->next;
1446 n_classes = 0;
1447 }
1448
1449 if(tree && tree->mod->ras->invar == NO) n_classes++;
1450
1451 if((tree->next && tree->next->is_mixt_tree == YES) || (!tree->next)) /*! current tree is the last element of this mixture */
1452 {
1453 if(n_classes < tree->mixt_tree->mod->ras->n_catg)
1454 {
1455 PhyML_Fprintf(stderr,"\n. %d class%s found in 'partitionelem' for file '%s' while",
1456 n_classes,
1457 (n_classes>1)?"es\0":"\0",
1458 tree->mixt_tree->io->in_align_file);
1459 PhyML_Fprintf(stderr,"\n. the corresponding 'siterates' element defined %d class%s.",
1460 tree->mixt_tree->mod->ras->n_catg,
1461 (tree->mixt_tree->mod->ras->n_catg>1)?"es\0":"\0");
1462 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
1463 Warn_And_Exit("\n");
1464 }
1465 }
1466
1467 tree = tree->next;
1468 }
1469 while(tree);
1470 }
1471 }
1472
1473 //////////////////////////////////////////////////////////////
1474 //////////////////////////////////////////////////////////////
1475
MIXT_Prune_Subtree(t_node * mixt_a,t_node * mixt_d,t_edge ** mixt_target,t_edge ** mixt_residual,t_tree * mixt_tree)1476 void MIXT_Prune_Subtree(t_node *mixt_a, t_node *mixt_d, t_edge **mixt_target, t_edge **mixt_residual, t_tree *mixt_tree)
1477 {
1478 t_node *a,*d;
1479 t_edge *target, *residual;
1480 t_tree *tree;
1481
1482 MIXT_Turn_Branches_OnOff_In_One_Elem(OFF,mixt_tree);
1483
1484 tree = mixt_tree;
1485 a = mixt_a;
1486 d = mixt_d;
1487 target = (mixt_target) ? (*mixt_target) : NULL;
1488 residual = (mixt_residual) ? (*mixt_residual) : NULL;
1489
1490 do
1491 {
1492 if(tree->is_mixt_tree == YES)
1493 {
1494 tree = tree->next;
1495 a = a->next;
1496 d = d->next;
1497 target = target ? target->next : NULL;
1498 residual = residual ? residual->next : NULL;
1499 }
1500
1501 Prune_Subtree(a,d,&target,&residual,tree);
1502
1503 tree = tree->next;
1504 a = a->next;
1505 d = d->next;
1506 target = target ? target->next : NULL;
1507 residual = residual ? residual->next : NULL;
1508 }
1509 while(tree && tree->is_mixt_tree == NO);
1510
1511 if(tree) Prune_Subtree(a,d,&target,&residual,tree);
1512
1513 /*! Turn branches of this mixt_tree to ON after recursive call
1514 to Prune_Subtree such that, if branches of mixt_tree->next
1515 point to those of mixt_tree, they are set to OFF when calling
1516 Prune */
1517 MIXT_Turn_Branches_OnOff_In_One_Elem(ON,mixt_tree);
1518
1519 }
1520
1521 //////////////////////////////////////////////////////////////
1522 //////////////////////////////////////////////////////////////
1523
MIXT_Graft_Subtree(t_edge * mixt_target,t_node * mixt_link,t_node * mixt_link_daughter,t_edge * mixt_residual,t_node * mixt_target_nd,t_tree * mixt_tree)1524 void MIXT_Graft_Subtree(t_edge *mixt_target, t_node *mixt_link, t_node *mixt_link_daughter, t_edge *mixt_residual, t_node *mixt_target_nd, t_tree *mixt_tree)
1525 {
1526 t_edge *target,*residual;
1527 t_node *link,*link_daughter,*target_nd;
1528 t_tree *tree;
1529
1530 MIXT_Turn_Branches_OnOff_In_One_Elem(OFF,mixt_tree);
1531
1532 tree = mixt_tree;
1533 target = mixt_target;
1534 residual = mixt_residual;
1535 link = mixt_link;
1536 link_daughter = mixt_link_daughter;
1537 target_nd = mixt_target_nd;
1538
1539 do
1540 {
1541 if(tree->is_mixt_tree == YES)
1542 {
1543 tree = tree->next;
1544 target = target->next;
1545 residual = residual->next;
1546 link = link->next;
1547 link_daughter = link_daughter ? link_daughter->next : NULL;
1548 target_nd = target_nd ? target_nd->next : NULL;
1549 }
1550
1551 Graft_Subtree(target,link,link_daughter,residual,target_nd,tree);
1552
1553 tree = tree->next;
1554 target = target->next;
1555 residual = residual->next;
1556 link = link->next;
1557 link_daughter = link_daughter ? link_daughter->next : NULL;
1558 target_nd = target_nd ? target_nd->next : NULL;
1559 }
1560 while(tree && tree->is_mixt_tree == NO);
1561
1562 if(tree) Graft_Subtree(target,link,link_daughter,residual,target_nd,tree);
1563
1564 /*! Turn branches of this mixt_tree to ON after recursive call
1565 to Graft_Subtree such that, if branches of mixt_tree->next
1566 point to those of mixt_tree, they are set to OFF when calling
1567 Graft */
1568 MIXT_Turn_Branches_OnOff_In_One_Elem(ON,mixt_tree);
1569 }
1570
1571 //////////////////////////////////////////////////////////////
1572 //////////////////////////////////////////////////////////////
1573
MIXT_Multiply_Scalar_Dbl(scalar_dbl * this,phydbl scalar)1574 void MIXT_Multiply_Scalar_Dbl(scalar_dbl *this, phydbl scalar)
1575 {
1576 if(this == NULL) return;
1577 else
1578 {
1579 scalar_dbl *buff;
1580 buff = this;
1581 do
1582 {
1583 buff->v *= scalar;
1584 buff = buff->next;
1585 }
1586 while(buff);
1587 }
1588 }
1589
1590 //////////////////////////////////////////////////////////////
1591 //////////////////////////////////////////////////////////////
1592
MIXT_Br_Len_Opt(t_edge * mixt_b,t_tree * mixt_tree)1593 void MIXT_Br_Len_Opt(t_edge *mixt_b, t_tree *mixt_tree)
1594 {
1595 t_edge *b;
1596 t_tree *tree;
1597 scalar_dbl *l;
1598
1599 MIXT_Turn_Branches_OnOff_In_All_Elem(ON,mixt_tree);
1600
1601 mixt_tree->ignore_mixt_info = YES;
1602
1603 b = mixt_b;
1604 tree = mixt_tree;
1605 l = NULL;
1606
1607 do
1608 {
1609 while(tree->is_mixt_tree == YES)
1610 {
1611 tree = tree->next;
1612 b = b->next;
1613 }
1614
1615 l = b->l;
1616 do
1617 {
1618 if(l->onoff == ON)
1619 {
1620 Br_Len_Opt(&(l->v),mixt_b,mixt_tree);
1621 l->onoff = OFF;
1622 }
1623 l = l->next;
1624 }
1625 while(l);
1626
1627 tree = tree->next;
1628 b = b->next;
1629 }
1630 while(tree);
1631
1632 mixt_tree->ignore_mixt_info = NO;
1633 }
1634
1635
1636 /*////////////////////////////////////////////////////////////
1637 ////////////////////////////////////////////////////////////*/
1638
MIXT_Br_Len_Involving_Invar(t_tree * mixt_tree)1639 void MIXT_Br_Len_Involving_Invar(t_tree *mixt_tree)
1640 {
1641 int i;
1642 scalar_dbl *l;
1643
1644 for(i=0;i<2*mixt_tree->n_otu-1;++i)
1645 {
1646 l = mixt_tree->a_edges[i]->l;
1647 do
1648 {
1649 l->v *= (1.-mixt_tree->mod->ras->pinvar->v);
1650 l = l->next;
1651 }
1652 while(l);
1653 }
1654 }
1655
1656 //////////////////////////////////////////////////////////////
1657 //////////////////////////////////////////////////////////////
1658
MIXT_Br_Len_Not_Involving_Invar(t_tree * mixt_tree)1659 void MIXT_Br_Len_Not_Involving_Invar(t_tree *mixt_tree)
1660 {
1661 int i;
1662 scalar_dbl *l;
1663
1664 for(i=0;i<2*mixt_tree->n_otu-1;++i)
1665 {
1666 l = mixt_tree->a_edges[i]->l;
1667 do
1668 {
1669 l->v /= (1.-mixt_tree->mod->ras->pinvar->v);
1670 l = l->next;
1671 }
1672 while(l);
1673 }
1674 }
1675
1676 //////////////////////////////////////////////////////////////
1677 //////////////////////////////////////////////////////////////
1678
MIXT_Unscale_Br_Len_Multiplier_Tree(t_tree * mixt_tree)1679 phydbl MIXT_Unscale_Br_Len_Multiplier_Tree(t_tree *mixt_tree)
1680 {
1681 int i;
1682 scalar_dbl *l;
1683
1684 for(i=0;i<2*mixt_tree->n_otu-1;++i)
1685 {
1686 l = mixt_tree->a_edges[i]->l;
1687 do
1688 {
1689 l->v /= mixt_tree->mod->br_len_mult->v;
1690 l = l->next;
1691 }
1692 while(l);
1693 }
1694 return(-1);
1695 }
1696
1697 //////////////////////////////////////////////////////////////
1698 //////////////////////////////////////////////////////////////
1699
MIXT_Rescale_Br_Len_Multiplier_Tree(t_tree * mixt_tree)1700 phydbl MIXT_Rescale_Br_Len_Multiplier_Tree(t_tree *mixt_tree)
1701 {
1702 int i;
1703 scalar_dbl *l;
1704
1705 for(i=0;i<2*mixt_tree->n_otu-1;++i)
1706 {
1707 l = mixt_tree->a_edges[i]->l;
1708 do
1709 {
1710 l->v *= mixt_tree->mod->br_len_mult->v;
1711 l = l->next;
1712 }
1713 while(l);
1714 }
1715 return(-1);
1716 }
1717
1718 //////////////////////////////////////////////////////////////
1719 //////////////////////////////////////////////////////////////
1720
MIXT_Rescale_Free_Rate_Tree(t_tree * mixt_tree)1721 phydbl MIXT_Rescale_Free_Rate_Tree(t_tree *mixt_tree)
1722 {
1723 int i,side_effect,at_boundary;
1724 t_edge *b;
1725
1726 side_effect = NO;
1727 for(i=0;i<2*mixt_tree->n_otu-1;++i)
1728 {
1729 b = mixt_tree->a_edges[i]->next;
1730
1731 at_boundary = NO;
1732 if(b->l->v > mixt_tree->mod->l_max-1.E-100 && b->l->v < mixt_tree->mod->l_max+1.E-100) at_boundary = YES;
1733 if(b->l->v > mixt_tree->mod->l_min-1.E-100 && b->l->v < mixt_tree->mod->l_min+1.E-100) at_boundary = YES;
1734
1735 b->l->v *= mixt_tree->mod->ras->free_rate_mr->v;
1736
1737 if(b->l->v > mixt_tree->mod->l_max && at_boundary == NO) side_effect = YES;
1738 if(b->l->v < mixt_tree->mod->l_min && at_boundary == NO) side_effect = YES;
1739 }
1740
1741 return side_effect;
1742 }
1743
1744 //////////////////////////////////////////////////////////////
1745 //////////////////////////////////////////////////////////////
1746
MIXT_Set_Ignore_Root(int yesno,t_tree * mixt_tree)1747 void MIXT_Set_Ignore_Root(int yesno, t_tree *mixt_tree)
1748 {
1749 t_tree *tree;
1750
1751 tree = mixt_tree;
1752 do
1753 {
1754 tree->ignore_root = yesno;
1755 tree = tree->next;
1756 }
1757 while(tree);
1758 }
1759
1760 //////////////////////////////////////////////////////////////
1761 //////////////////////////////////////////////////////////////
1762
MIXT_Set_Alias_Subpatt(int onoff,t_tree * mixt_tree)1763 void MIXT_Set_Alias_Subpatt(int onoff, t_tree *mixt_tree)
1764 {
1765 t_tree *tree;
1766
1767 tree = mixt_tree;
1768 do
1769 {
1770 tree->update_alias_subpatt = onoff;
1771 tree = tree->next;
1772 }
1773 while(tree);
1774 }
1775
1776 //////////////////////////////////////////////////////////////
1777 //////////////////////////////////////////////////////////////
1778
MIXT_Check_Edge_Lens_In_All_Elem(t_tree * mixt_tree)1779 void MIXT_Check_Edge_Lens_In_All_Elem(t_tree *mixt_tree)
1780 {
1781 t_tree *tree;
1782
1783 /*! Check that all the edges in a mixt_tree at pointing
1784 to a single set of lengths
1785 */
1786 tree = mixt_tree;
1787 do
1788 {
1789 MIXT_Check_Edge_Lens_In_One_Elem(tree);
1790 tree = tree->next_mixt;
1791 }
1792 while(tree);
1793 }
1794
1795 //////////////////////////////////////////////////////////////
1796 //////////////////////////////////////////////////////////////
1797
MIXT_Check_Edge_Lens_In_One_Elem(t_tree * mixt_tree)1798 void MIXT_Check_Edge_Lens_In_One_Elem(t_tree *mixt_tree)
1799 {
1800 t_tree *tree;
1801 int i;
1802
1803 tree = mixt_tree->next;
1804 do
1805 {
1806 if(tree->next && tree->next->is_mixt_tree == NO)
1807 {
1808 for(i=0;i<2*tree->n_otu-1;++i)
1809 {
1810 if(tree->a_edges[i]->l != tree->next->a_edges[i]->l)
1811 {
1812 PhyML_Fprintf(stderr,"\n. %p %p",tree->a_edges[i]->l,tree->next->a_edges[i]->l);
1813 PhyML_Fprintf(stderr,"\n. Only one set of edge lengths is allowed ");
1814 PhyML_Fprintf(stderr,"\n. in a 'partitionelem'. Please fix your XML file.");
1815 Exit("\n");
1816 }
1817 }
1818 }
1819 tree = tree->next;
1820 }
1821 while(tree && tree->next && tree->next->is_mixt_tree == NO);
1822 }
1823
1824 //////////////////////////////////////////////////////////////
1825 //////////////////////////////////////////////////////////////
1826
MIXT_Pars(t_edge * mixt_b,t_tree * mixt_tree)1827 int MIXT_Pars(t_edge *mixt_b, t_tree *mixt_tree)
1828 {
1829 t_edge *b;
1830 t_tree *tree;
1831
1832 b = mixt_b;
1833 tree = mixt_tree;
1834 mixt_tree->c_pars = 0;
1835
1836 do
1837 {
1838 if(tree->next)
1839 {
1840 Pars(b?b->next:NULL,tree->next);
1841 mixt_tree->c_pars += tree->next->c_pars;
1842 }
1843
1844 if(mixt_b != NULL) b = b->next_mixt;
1845 tree = tree->next_mixt;
1846 }
1847 while(tree);
1848
1849 tree = mixt_tree;
1850 do
1851 {
1852 tree->c_pars = mixt_tree->c_pars;
1853 tree = tree->next_mixt;
1854 }
1855 while(tree);
1856
1857 return(mixt_tree->c_pars);
1858 }
1859
1860 //////////////////////////////////////////////////////////////
1861 //////////////////////////////////////////////////////////////
1862
MIXT_Bootstrap(char * best_tree,xml_node * root)1863 void MIXT_Bootstrap(char *best_tree, xml_node *root)
1864 {
1865 xml_node *n,*p_elem;
1866 char *bootstrap;
1867
1868 assert(root);
1869
1870 n = XML_Search_Node_Name("phyml",NO,root);
1871
1872 bootstrap = XML_Get_Attribute_Value(n,"bootstrap");
1873
1874 if(!bootstrap) return;
1875 else
1876 {
1877 int n_boot,i,j,k;
1878 xml_attr *boot_attr,*seqfile_attr,*out_attr,*boot_out_attr;
1879 char *orig_align,*boot_out_file_name,*xml_boot_file_name,*buff;
1880 FILE *boot_fp_in_align,*xml_boot_file_fp;
1881 option *io,*dum;
1882 align **boot_data,**orig_data;
1883 int position,elem;
1884 xml_node *boot_root;
1885 int pid;
1886 char *s;
1887
1888 orig_align = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1889
1890 xml_boot_file_name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1891 strcpy(xml_boot_file_name,"phyml_boot_config.");
1892 pid = (int)getpid();
1893 sprintf(xml_boot_file_name+strlen(xml_boot_file_name),"%d",pid);
1894 strcat(xml_boot_file_name,".xml");
1895
1896 out_attr = XML_Search_Attribute(root,"output.file");
1897 assert(out_attr);
1898 boot_out_file_name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1899 strcpy(boot_out_file_name,out_attr->value);
1900 s = XML_Get_Attribute_Value(root,"run.id");
1901 if(s)
1902 {
1903 strcat(boot_out_file_name,"_");
1904 strcat(boot_out_file_name,s);
1905 }
1906
1907 n_boot = atoi(bootstrap);
1908
1909 io = NULL;
1910 for(i=0;i<n_boot;++i)
1911 {
1912 boot_root = XML_Copy_XML_Graph(root);
1913
1914 /*! Set the number of bootstrap repeats to 0
1915 in each generated XML file */
1916 boot_attr = XML_Search_Attribute(boot_root,"bootstrap");
1917 assert(boot_attr);
1918 strcpy(boot_attr->value,"0");
1919
1920 /*! Set the output file name for each bootstrap analysis */
1921 boot_out_attr = XML_Search_Attribute(boot_root,"output.file");
1922 assert(boot_out_attr);
1923 buff = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1924 strcpy(buff,boot_out_attr->value);
1925 Free(boot_out_attr->value);
1926 boot_out_attr->value = buff;
1927 sprintf(boot_out_attr->value+strlen(boot_out_attr->value),"_boot.%d",pid);
1928
1929 p_elem = boot_root;
1930 elem = 0;
1931 do
1932 {
1933 p_elem = XML_Search_Node_Name("partitionelem",YES,p_elem);
1934 if(!p_elem) break;
1935
1936 io = (option *)Make_Input();
1937 Set_Defaults_Input(io);
1938
1939 /*! Get the original sequence file name and the corresponding
1940 attribute in the XML graph
1941 */
1942 seqfile_attr = NULL;
1943 seqfile_attr = XML_Search_Attribute(p_elem,"file.name");
1944 assert(seqfile_attr);
1945
1946 strcpy(orig_align,seqfile_attr->value);
1947
1948 /*! Open the original sequence file */
1949 io->fp_in_align = Openfile(orig_align,0);
1950
1951 /*! Read in the original sequence file */
1952 orig_data = Get_Seq(io);
1953 rewind(io->fp_in_align);
1954
1955 /*! Read in the original sequence file and put
1956 it in 'boot_data' structure */
1957 boot_data = Get_Seq(io);
1958
1959 fclose(io->fp_in_align);
1960
1961 /*! Bootstrap resampling: sample from original and put in boot */
1962 for(j=0;j<boot_data[0]->len;++j)
1963 {
1964 position = Rand_Int(0,(int)(boot_data[0]->len-1.0));
1965 for(k=0;k<io->n_otu;++k) boot_data[k]->state[j] = orig_data[k]->state[position];
1966 }
1967
1968 /*! Modify the sequence file attribute in the original XML
1969 graph */
1970 buff = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1971 Free(seqfile_attr->value);
1972 seqfile_attr->value = buff;
1973 sprintf(seqfile_attr->value,"%s_%d_%d",orig_align,elem,i);
1974
1975 /*! Open a new sequence file with the modified attribute name */
1976 boot_fp_in_align = Openfile(seqfile_attr->value,1);
1977
1978 /*! Print the bootstrap data set in it */
1979 Print_Seq(boot_fp_in_align,boot_data,io->n_otu);
1980 fclose(boot_fp_in_align);
1981
1982 Free_Seq(orig_data,io->n_otu);
1983 Free_Seq(boot_data,io->n_otu);
1984
1985 Free_Input(io);
1986 elem++;
1987 }
1988 while(p_elem);
1989
1990 /*! Open bootstrap XML file in writing mode */
1991 xml_boot_file_fp = Openfile(xml_boot_file_name,1);
1992
1993 /*! Write the bootstrap XML graph */
1994 XML_Write_XML_Graph(xml_boot_file_fp,boot_root);
1995 fclose(xml_boot_file_fp);
1996
1997 /*! Reconstruct the tree */
1998 dum = PhyML_XML(xml_boot_file_name);
1999 Free(dum);
2000
2001 /*! Remove the bootstrap alignment files */
2002 p_elem = boot_root;
2003 do
2004 {
2005 p_elem = XML_Search_Node_Name("partitionelem",YES,p_elem);
2006 if(!p_elem) break;
2007 seqfile_attr = XML_Search_Attribute(p_elem,"file.name");
2008 unlink(seqfile_attr->value);
2009 }
2010 while(p_elem);
2011
2012 XML_Free_XML_Tree(boot_root);
2013 }
2014
2015 Free(xml_boot_file_name);
2016 Free(orig_align);
2017 Free(boot_out_file_name);
2018 }
2019
2020 }
2021
2022 //////////////////////////////////////////////////////////////
2023 //////////////////////////////////////////////////////////////
2024
MIXT_Set_Pars_Thresh(t_tree * mixt_tree)2025 void MIXT_Set_Pars_Thresh(t_tree *mixt_tree)
2026 {
2027 t_tree *tree;
2028
2029 tree = mixt_tree;
2030 do
2031 {
2032 tree->mod->s_opt->pars_thresh = (tree->io->datatype == AA)?(15):(5);
2033 tree = tree->next_mixt;
2034 }
2035 while(tree);
2036 }
2037
2038 //////////////////////////////////////////////////////////////
2039 //////////////////////////////////////////////////////////////
2040
MIXT_Get_Mean_Edge_Len(t_edge * mixt_b,t_tree * mixt_tree)2041 phydbl MIXT_Get_Mean_Edge_Len(t_edge *mixt_b, t_tree *mixt_tree)
2042 {
2043 phydbl sum;
2044 int n;
2045 t_tree *tree;
2046 t_edge *b;
2047
2048 if(mixt_tree->is_mixt_tree == NO) return mixt_b->l->v;
2049
2050 b = mixt_b;
2051 tree = mixt_tree;
2052 sum = .0;
2053 n = 0 ;
2054 do
2055 {
2056 if(tree->is_mixt_tree == YES)
2057 {
2058 tree = tree->next;
2059 b = b->next;
2060 }
2061
2062 sum += b->l->v * (tree->mixt_tree ? tree->mixt_tree->mod->br_len_mult->v : 1.0);
2063 n++;
2064 b = b->next;
2065 tree = tree->next;
2066 }
2067 while(b);
2068
2069 return(sum / (phydbl)n);
2070 }
2071
2072 //////////////////////////////////////////////////////////////
2073 //////////////////////////////////////////////////////////////
2074
MIXT_Get_Sum_Chained_Scalar_Dbl(scalar_dbl * s)2075 phydbl MIXT_Get_Sum_Chained_Scalar_Dbl(scalar_dbl *s)
2076 {
2077 scalar_dbl *s_buff;
2078 phydbl sum;
2079
2080 s_buff = s;
2081 sum = .0;
2082 do
2083 {
2084 sum += s_buff->v;
2085 s_buff = s_buff->next;
2086 }
2087 while(s_buff);
2088
2089 return sum;
2090
2091 }
2092
2093 //////////////////////////////////////////////////////////////
2094 //////////////////////////////////////////////////////////////
2095
MIXT_Get_Sum_Of_Probas_Across_Mixtures(phydbl r_mat_weight_sum,phydbl e_frq_weight_sum,t_tree * mixt_tree)2096 phydbl MIXT_Get_Sum_Of_Probas_Across_Mixtures(phydbl r_mat_weight_sum, phydbl e_frq_weight_sum, t_tree *mixt_tree)
2097 {
2098 t_tree *tree;
2099 phydbl sum;
2100
2101 sum = .0;
2102 tree = mixt_tree->next;
2103 do
2104 {
2105 // e.g., if mixture has two classes, one of these
2106 // corresponding to invariable sites. We need to skip it.
2107 if(tree->mod->ras->invar == YES) tree = tree->next;
2108
2109 sum +=
2110 mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number] *
2111 tree->mod->r_mat_weight->v / r_mat_weight_sum *
2112 tree->mod->e_frq_weight->v / e_frq_weight_sum;
2113
2114
2115 tree = tree->next;
2116
2117 }
2118 while(tree && tree->is_mixt_tree == NO);
2119
2120 return(sum);
2121 }
2122
2123 //////////////////////////////////////////////////////////////
2124 //////////////////////////////////////////////////////////////
2125
MIXT_Set_Br_Len_Var(t_edge * mixt_b,t_tree * mixt_tree)2126 void MIXT_Set_Br_Len_Var(t_edge *mixt_b, t_tree *mixt_tree)
2127 {
2128 t_tree *tree;
2129
2130 if(mixt_b != NULL)
2131 {
2132 t_edge *b;
2133
2134 tree = mixt_tree->next;
2135 b = mixt_b->next;
2136
2137 do
2138 {
2139 Set_Br_Len_Var(b,tree);
2140 tree = tree->next;
2141 b = b->next;
2142 }
2143 while(tree);
2144 }
2145 else
2146 {
2147 tree = mixt_tree->next;
2148
2149 do
2150 {
2151 Set_Br_Len_Var(NULL,tree);
2152 tree = tree->next;
2153 }
2154 while(tree);
2155
2156 }
2157 }
2158
2159 //////////////////////////////////////////////////////////////
2160 //////////////////////////////////////////////////////////////
2161
MIXT_Set_Br_Len(phydbl val,t_edge * mixt_b,t_tree * mixt_tree)2162 void MIXT_Set_Br_Len(phydbl val, t_edge *mixt_b, t_tree *mixt_tree)
2163 {
2164 scalar_dbl *l;
2165
2166 l = mixt_b->l;
2167 do
2168 {
2169 l->v = val;
2170 l = l->next;
2171 }
2172 while(l);
2173 }
2174
2175 //////////////////////////////////////////////////////////////
2176 //////////////////////////////////////////////////////////////
2177
MIXT_Add_Root(t_edge * mixt_b,t_tree * mixt_tree)2178 void MIXT_Add_Root(t_edge *mixt_b, t_tree *mixt_tree)
2179 {
2180 t_tree *tree;
2181 t_edge *b;
2182
2183 tree = mixt_tree;
2184 b = mixt_b;
2185 do
2186 {
2187 if(tree->is_mixt_tree)
2188 {
2189 tree = tree->next;
2190 b = b->next;
2191 }
2192
2193 // Condition is true when tree is not chained
2194 if(b == NULL) break;
2195
2196 Add_Root(b,tree);
2197
2198 tree = tree->next;
2199 b = b->next;
2200 }
2201 while(tree);
2202
2203 tree = mixt_tree;
2204 do
2205 {
2206 assert(tree->n_root != NULL);
2207
2208 if(tree->next) tree->n_root->next = tree->next->n_root;
2209 if(tree->prev) tree->n_root->prev = tree->prev->n_root;
2210 if(tree->next_mixt) tree->n_root->next_mixt = tree->next_mixt->n_root;
2211 if(tree->prev_mixt) tree->n_root->prev_mixt = tree->prev_mixt->n_root;
2212
2213 tree = tree->next;
2214 }
2215 while(tree);
2216
2217 }
2218
2219 //////////////////////////////////////////////////////////////
2220 //////////////////////////////////////////////////////////////
2221
MIXT_RATES_Update_Cur_Bl(t_tree * mixt_tree)2222 void MIXT_RATES_Update_Cur_Bl(t_tree *mixt_tree)
2223 {
2224 t_tree *tree;
2225
2226 tree = mixt_tree->next;
2227 do
2228 {
2229 RATES_Update_Cur_Bl(tree);
2230 tree = tree->next;
2231 }
2232 while(tree);
2233 }
2234
2235 //////////////////////////////////////////////////////////////
2236 //////////////////////////////////////////////////////////////
2237
MIXT_Update_Br_Len_Multipliers(t_mod * mod)2238 void MIXT_Update_Br_Len_Multipliers(t_mod *mod)
2239 {
2240 phydbl sum;
2241 t_mod *loc;
2242 int n_mixt;
2243
2244 loc = mod;
2245 sum = 0.0;
2246 n_mixt = 0;
2247 do
2248 {
2249 /* if(loc->s_opt->opt_br_len_mult == YES) */
2250 /* { */
2251 sum += loc->br_len_mult_unscaled->v;
2252 n_mixt++;
2253 /* } */
2254 loc = loc->next_mixt;
2255 }
2256 while(loc);
2257
2258 loc = mod;
2259 do
2260 {
2261 if(loc->s_opt->opt_br_len_mult == YES)
2262 {
2263 loc->br_len_mult->v = loc->br_len_mult_unscaled->v / sum;
2264 loc->br_len_mult->v *= (phydbl)(n_mixt);
2265 /* printf("\n. HERE %f %f\n",loc->br_len_mult_unscaled->v,loc->br_len_mult->v); */
2266 }
2267 loc = loc->next_mixt;
2268 }
2269 while(loc);
2270 }
2271
2272 //////////////////////////////////////////////////////////////
2273 //////////////////////////////////////////////////////////////
2274
MIXT_Init_Model(t_tree * mixt_tree)2275 void MIXT_Init_Model(t_tree *mixt_tree)
2276 {
2277 t_mod *mod,*mixt_mod;
2278 option *io;
2279 t_tree *tree;
2280
2281 assert(mixt_tree);
2282
2283 mixt_mod = mixt_tree->mod;
2284 io = mixt_tree->io;
2285
2286 mod = mixt_mod;
2287 do
2288 {
2289 Init_Model(mod->io->cdata,mod,io);
2290 mod = mod->next;
2291 }
2292 while(mod);
2293
2294 tree = mixt_tree;
2295 do
2296 {
2297 if(tree->next_mixt != NULL)
2298 {
2299 tree->mod->next_mixt = tree->next_mixt->mod;
2300 tree->next_mixt->mod->prev_mixt = tree->mod;
2301 }
2302 tree = tree->next_mixt;
2303 }
2304 while(tree);
2305 }
2306
2307 //////////////////////////////////////////////////////////////
2308 //////////////////////////////////////////////////////////////
2309
MIXT_Check_Model_Validity(t_tree * mixt_tree)2310 void MIXT_Check_Model_Validity(t_tree *mixt_tree)
2311 {
2312 // Verify that models associated to distinct data partition elements do not
2313 // share the same empirical character frequencies.
2314
2315 t_mod *mod_in, *mod_out;
2316
2317 mod_out = mixt_tree->mod;
2318
2319 do
2320 {
2321 mod_in = mod_out;
2322 do
2323 {
2324 if(mod_in->io->cdata != mod_out->io->cdata)
2325 {
2326 if(mod_in->e_frq == mod_out->e_frq)
2327 {
2328 if(mod_in->io->datatype == NT &&
2329 mod_in->e_frq->user_state_freq == NO &&
2330 mod_in->whichmodel != JC69 &&
2331 mod_in->whichmodel != K80)
2332 {
2333 PhyML_Fprintf(stderr,"\n. A vector of observed nucleotide frequencies should correspond ");
2334 PhyML_Fprintf(stderr,"\n. to one data set only. If you are using the XML interface, ");
2335 PhyML_Fprintf(stderr,"\n. please amend your file accordingly.");
2336 Exit("\n");
2337 }
2338 else if(mod_in->io->datatype == AA && mod_in->e_frq->empirical_state_freq == YES)
2339 {
2340 PhyML_Fprintf(stderr,"\n. A vector of observed amino-acid frequencies should correspond ");
2341 PhyML_Fprintf(stderr,"\n. to one data set only. If you are using the XML interface, ");
2342 PhyML_Fprintf(stderr,"\n. please amend your file accordingly.");
2343 Exit("\n");
2344 }
2345 }
2346 }
2347 mod_in = mod_in->next;
2348 }
2349 while(mod_in);
2350 mod_out = mod_out->next;
2351 }
2352 while(mod_out);
2353
2354
2355 }
2356
2357 //////////////////////////////////////////////////////////////
2358 //////////////////////////////////////////////////////////////
2359
MIXT_Starting_Tree(t_tree * mixt_tree)2360 t_tree *MIXT_Starting_Tree(t_tree *mixt_tree)
2361 {
2362 t_tree *tree;
2363
2364 tree = NULL;
2365
2366 if(mixt_tree->io->mod->s_opt->random_input_tree == NO)
2367 {
2368 switch(mixt_tree->io->in_tree)
2369 {
2370 case 2: // user-defined input tree
2371 {
2372 assert(mixt_tree->io->fp_in_tree);
2373
2374 tree = Read_User_Tree(mixt_tree->io->cdata,
2375 mixt_tree->mod,
2376 mixt_tree->io);
2377
2378 break;
2379 }
2380 case 1: case 0:
2381 {
2382 // Build a BioNJ tree from the analysis of
2383 // the first partition element
2384 tree = Dist_And_BioNJ(mixt_tree->data,
2385 mixt_tree->mod,
2386 mixt_tree->io);
2387 break;
2388 }
2389 default : assert(FALSE);
2390 }
2391 }
2392 else
2393 {
2394 tree = Make_Tree_From_Scratch(mixt_tree->n_otu,mixt_tree->data);
2395 Random_Tree(tree);
2396 }
2397
2398 return tree;
2399 }
2400
2401 //////////////////////////////////////////////////////////////
2402 //////////////////////////////////////////////////////////////
2403
MIXT_Connect_Cseqs_To_Nodes(t_tree * mixt_tree)2404 void MIXT_Connect_Cseqs_To_Nodes(t_tree *mixt_tree)
2405 {
2406 t_tree *tree;
2407
2408 Copy_Tree(mixt_tree,mixt_tree->next);
2409
2410 tree = mixt_tree;
2411 do
2412 {
2413 Connect_CSeqs_To_Nodes(tree->data,mixt_tree->io,tree);
2414 tree = tree->next;
2415 }
2416 while(tree);
2417
2418 }
2419
2420 //////////////////////////////////////////////////////////////
2421 //////////////////////////////////////////////////////////////
2422
MIXT_Init_T_Beg(t_tree * mixt_tree)2423 void MIXT_Init_T_Beg(t_tree *mixt_tree)
2424 {
2425 t_tree *tree;
2426
2427 /*! Initialize t_beg in each mixture tree */
2428 tree = mixt_tree;
2429 do
2430 {
2431 time(&(tree->t_beg));
2432 tree = tree->next_mixt;
2433 }
2434 while(tree);
2435 }
2436
2437 //////////////////////////////////////////////////////////////
2438 //////////////////////////////////////////////////////////////
2439
MIXT_Init_T_End(t_tree * mixt_tree)2440 void MIXT_Init_T_End(t_tree *mixt_tree)
2441 {
2442 t_tree *tree;
2443
2444 /*! Initialize t_beg in each mixture tree */
2445 tree = mixt_tree;
2446 do
2447 {
2448 time(&(tree->t_current));
2449 tree = tree->next_mixt;
2450 }
2451 while(tree);
2452 }
2453
2454 //////////////////////////////////////////////////////////////
2455 //////////////////////////////////////////////////////////////
2456
MIXT_Prepare_All(int num_rand_tree,t_tree * mixt_tree)2457 void MIXT_Prepare_All(int num_rand_tree, t_tree *mixt_tree)
2458 {
2459 t_tree *tree;
2460
2461 MIXT_Check_Model_Validity(mixt_tree);
2462 MIXT_Init_Model(mixt_tree);
2463 Print_Data_Structure(NO,stdout,mixt_tree);
2464 tree = MIXT_Starting_Tree(mixt_tree);
2465 Copy_Tree(tree,mixt_tree);
2466 Free_Tree(tree);
2467
2468 if(mixt_tree->io->mod->s_opt->random_input_tree)
2469 {
2470 PhyML_Printf("\n\n. [%3d/%3d]",num_rand_tree+1,mixt_tree->io->mod->s_opt->n_rand_starts);
2471 Random_Tree(mixt_tree);
2472 }
2473
2474 MIXT_Connect_Cseqs_To_Nodes(mixt_tree);
2475 MIXT_Init_T_Beg(mixt_tree);
2476
2477 MIXT_Make_Tree_For_Pars(mixt_tree);
2478 MIXT_Make_Tree_For_Lk(mixt_tree);
2479 MIXT_Make_Spr(mixt_tree);
2480
2481 MIXT_Chain_All(mixt_tree);
2482 MIXT_Check_Edge_Lens_In_All_Elem(mixt_tree);
2483 MIXT_Turn_Branches_OnOff_In_All_Elem(ON,mixt_tree);
2484 MIXT_Check_Invar_Struct_In_Each_Partition_Elem(mixt_tree);
2485 MIXT_Check_RAS_Struct_In_Each_Partition_Elem(mixt_tree);
2486 }
2487
2488 //////////////////////////////////////////////////////////////
2489 //////////////////////////////////////////////////////////////
2490
MIXT_Make_Spr(t_tree * mixt_tree)2491 void MIXT_Make_Spr(t_tree *mixt_tree)
2492 {
2493 t_tree *tree;
2494
2495 tree = mixt_tree;
2496 do
2497 {
2498 Make_Spr(tree);
2499 tree = tree->next;
2500 }
2501 while(tree);
2502 }
2503
2504 //////////////////////////////////////////////////////////////
2505 //////////////////////////////////////////////////////////////
2506
MIXT_Make_Tree_For_Pars(t_tree * mixt_tree)2507 void MIXT_Make_Tree_For_Pars(t_tree *mixt_tree)
2508 {
2509 t_tree *tree;
2510
2511 tree = mixt_tree;
2512 do
2513 {
2514 if(tree->is_mixt_tree == NO) Make_Tree_For_Pars(tree);
2515 tree = tree->next;
2516 }
2517 while(tree);
2518 }
2519
2520 //////////////////////////////////////////////////////////////
2521 //////////////////////////////////////////////////////////////
2522
MIXT_Make_Tree_For_Lk(t_tree * mixt_tree)2523 void MIXT_Make_Tree_For_Lk(t_tree *mixt_tree)
2524 {
2525 t_tree *tree;
2526
2527 tree = mixt_tree;
2528 do
2529 {
2530 if(tree->is_mixt_tree == NO) Make_Tree_For_Lk(tree);
2531 else
2532 {
2533 tree->c_lnL_sorted = (phydbl *)mCalloc(tree->n_pattern,sizeof(phydbl));
2534 tree->cur_site_lk = (phydbl *)mCalloc(tree->n_pattern,sizeof(phydbl));
2535 tree->old_site_lk = (phydbl *)mCalloc(tree->n_pattern,sizeof(phydbl));
2536 tree->site_lk_cat = (phydbl *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(phydbl));
2537 tree->unscaled_site_lk_cat = (phydbl *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->n_pattern,sizeof(phydbl));
2538 tree->fact_sum_scale = (int *)mCalloc(tree->n_pattern,sizeof(int));
2539
2540
2541 #if (defined(__AVX__) || defined(__SSE3__))
2542 #ifndef WIN32
2543 if(posix_memalign((void **)&tree->expl,BYTE_ALIGN,(size_t)2*tree->mod->n_mixt_classes*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
2544 #else
2545 tree->expl = _aligned_malloc(2*tree->mod->n_mixt_classes*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
2546 #endif
2547 #else
2548 tree->expl = (phydbl *)mCalloc(2*tree->mod->n_mixt_classes*tree->mod->ns,sizeof(phydbl));
2549 #endif
2550
2551 for(int i=0;i<2*tree->n_otu-1;++i) Make_Edge_NNI(tree->a_edges[i]);
2552
2553 tree->log_lks_aLRT = (phydbl **)mCalloc(3,sizeof(phydbl *));
2554 for(int i=0;i<3;i++) tree->log_lks_aLRT[i] = (phydbl *)mCalloc(tree->data->init_len,sizeof(phydbl));
2555 }
2556
2557 tree = tree->next;
2558 }
2559 while(tree);
2560 }
2561
2562 //////////////////////////////////////////////////////////////
2563 //////////////////////////////////////////////////////////////
2564
MIXT_Ancestral_Sequences_One_Node(t_node * mixt_d,t_tree * mixt_tree,int print)2565 void MIXT_Ancestral_Sequences_One_Node(t_node *mixt_d, t_tree *mixt_tree, int print)
2566 {
2567 if(mixt_d->tax) return;
2568 else
2569 {
2570 t_node *v0,*v1,*v2; // three neighbours of d
2571 t_edge *b0,*b1,*b2;
2572 int i,j;
2573 int catg;
2574 phydbl p0, p1, p2;
2575 phydbl *p;
2576 t_node *d,*curr_mixt_d;
2577 t_tree *tree, *curr_mixt_tree;
2578 int site,csite;
2579 phydbl *p_lk0, *p_lk1, *p_lk2;
2580 int *sum_scale0, *sum_scale1, *sum_scale2;
2581 phydbl r_mat_weight_sum, e_frq_weight_sum, sum_probas;
2582 phydbl *Pij0, *Pij1, *Pij2;
2583 int NsNs, Ns, NsNg;
2584 FILE *fp;
2585
2586 if(!mixt_d) return;
2587
2588 curr_mixt_tree = mixt_tree;
2589 curr_mixt_d = mixt_d;
2590 fp = mixt_tree->io->fp_out_ancestral_seq;
2591
2592
2593 do /* For each partition element */
2594 {
2595 if(curr_mixt_tree->next)
2596 {
2597 r_mat_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(curr_mixt_tree->next->mod->r_mat_weight);
2598 e_frq_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(curr_mixt_tree->next->mod->e_frq_weight);
2599 sum_probas = MIXT_Get_Sum_Of_Probas_Across_Mixtures(r_mat_weight_sum, e_frq_weight_sum, curr_mixt_tree);
2600 }
2601 else
2602 {
2603 r_mat_weight_sum = 1.;
2604 e_frq_weight_sum = 1.;
2605 sum_probas = 1.;
2606 }
2607
2608 Ns = curr_mixt_tree->next ? curr_mixt_tree->next->mod->ns : curr_mixt_tree->mod->ns;
2609 NsNs = Ns*Ns;
2610 NsNg = Ns*curr_mixt_tree->mod->ras->n_catg;
2611
2612 p = (phydbl *)mCalloc(Ns,sizeof(phydbl));
2613
2614 /* for(site=0;site<curr_mixt_tree->n_pattern;site++) // For each site in the current partition element */
2615 for(site=0;site<curr_mixt_tree->data->init_len;site++) // For each site in the current partition element
2616 {
2617 csite = curr_mixt_tree->data->sitepatt[site];
2618
2619 d = curr_mixt_d->next ? curr_mixt_d->next : curr_mixt_d;
2620 tree = curr_mixt_tree->next ? curr_mixt_tree->next : curr_mixt_tree;
2621
2622 for(i=0;i<tree->mod->ns;i++) p[i] = .0;
2623
2624 do // For each class of the mixture model that applies to the current partition element
2625 {
2626 if(tree->is_mixt_tree == YES)
2627 {
2628 tree = tree->next;
2629 d = d->next;
2630 }
2631
2632 v0 = d->v[0];
2633 v1 = d->v[1];
2634 v2 = d->v[2];
2635
2636 b0 = d->b[0];
2637 b1 = d->b[1];
2638 b2 = d->b[2];
2639
2640 Pij0 = b0->Pij_rr;
2641 Pij1 = b1->Pij_rr;
2642 Pij2 = b2->Pij_rr;
2643
2644 if(v0 == b0->left)
2645 {
2646 p_lk0 = b0->p_lk_left;
2647 sum_scale0 = b0->sum_scale_left;
2648 }
2649 else
2650 {
2651 p_lk0 = b0->p_lk_rght;
2652 sum_scale0 = b0->sum_scale_rght;
2653 }
2654
2655 if(v1 == b1->left)
2656 {
2657 p_lk1 = b1->p_lk_left;
2658 sum_scale1 = b1->sum_scale_left;
2659 }
2660 else
2661 {
2662 p_lk1 = b1->p_lk_rght;
2663 sum_scale1 = b1->sum_scale_rght;
2664 }
2665
2666 if(v2 == b2->left)
2667 {
2668 p_lk2 = b2->p_lk_left;
2669 sum_scale2 = b2->sum_scale_left;
2670 }
2671 else
2672 {
2673 p_lk2 = b2->p_lk_rght;
2674 sum_scale2 = b2->sum_scale_rght;
2675 }
2676
2677
2678 for(catg=0;catg<tree->mod->ras->n_catg;catg++)
2679 {
2680 for(i=0;i<Ns;i++)
2681 {
2682 p0 = .0;
2683 if(v0->tax)
2684 for(j=0;j<tree->mod->ns;j++)
2685 {
2686 p0 += v0->b[0]->p_lk_tip_r[csite*Ns+j] * Pij0[catg*NsNs+i*Ns+j];
2687
2688 /* printf("\n. p0 %d %f", */
2689 /* v0->b[0]->p_lk_tip_r[site*Ns+j], */
2690 /* Pij0[catg*NsNs+i*Ns+j]); */
2691 }
2692 else
2693 for(j=0;j<tree->mod->ns;j++)
2694 {
2695 p0 += p_lk0[csite*NsNg+catg*Ns+j] * Pij0[catg*NsNs+i*Ns+j] / (phydbl)POW(2,sum_scale0[catg*curr_mixt_tree->n_pattern+csite]);
2696
2697 /* p0 += p_lk0[site*NsNg+catg*Ns+j] * Pij0[catg*NsNs+i*Ns+j]; */
2698
2699 /* printf("\n. p0 %f %f", */
2700 /* p_lk0[site*NsNg+catg*Ns+j], */
2701 /* Pij0[catg*NsNs+i*Ns+j]); */
2702 }
2703 p1 = .0;
2704 if(v1->tax)
2705 for(j=0;j<tree->mod->ns;j++)
2706 {
2707 p1 += v1->b[0]->p_lk_tip_r[csite*Ns+j] * Pij1[catg*NsNs+i*Ns+j];
2708
2709 /* printf("\n. p1 %d %f", */
2710 /* v1->b[0]->p_lk_tip_r[site*Ns+j], */
2711 /* Pij1[catg*NsNs+i*Ns+j]); */
2712 }
2713
2714 else
2715 for(j=0;j<tree->mod->ns;j++)
2716 {
2717 p1 += p_lk1[csite*NsNg+catg*Ns+j] * Pij1[catg*NsNs+i*Ns+j] / (phydbl)POW(2,sum_scale1[catg*curr_mixt_tree->n_pattern+csite]);
2718
2719 /* p1 += p_lk1[site*NsNg+catg*Ns+j] * Pij1[catg*NsNs+i*Ns+j]; */
2720
2721 /* printf("\n. p1 %f %f", */
2722 /* p_lk1[site*NsNg+catg*Ns+j], */
2723 /* Pij1[catg*NsNs+i*Ns+j]); */
2724 }
2725
2726
2727 p2 = .0;
2728 if(v2->tax)
2729 for(j=0;j<tree->mod->ns;j++)
2730 {
2731 p2 += v2->b[0]->p_lk_tip_r[csite*Ns+j] * Pij2[catg*NsNs+i*Ns+j];
2732 /* printf("\n. p2 %d %f", */
2733 /* v2->b[0]->p_lk_tip_r[site*Ns+j], */
2734 /* Pij2[catg*NsNs+i*Ns+j]); */
2735 }
2736 else
2737 for(j=0;j<tree->mod->ns;j++)
2738 {
2739 p2 += p_lk2[csite*NsNg+catg*Ns+j] * Pij2[catg*NsNs+i*Ns+j] / (phydbl)POW(2,sum_scale2[catg*curr_mixt_tree->n_pattern+csite]);
2740
2741 /* p2 += p_lk2[site*NsNg+catg*Ns+j] * Pij2[catg*NsNs+i*Ns+j]; */
2742
2743 /* printf("\n. p2 %f %f", */
2744 /* p_lk2[site*NsNg+catg*Ns+j], */
2745 /* Pij2[catg*NsNs+i*Ns+j]); */
2746 }
2747
2748 p[i] +=
2749 p0*p1*p2*
2750 tree->mod->e_frq->pi->v[i] /
2751 tree->cur_site_lk[csite] *
2752 curr_mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number] *
2753 tree->mod->r_mat_weight->v / r_mat_weight_sum *
2754 tree->mod->e_frq_weight->v / e_frq_weight_sum /
2755 sum_probas;
2756
2757 if(print == YES)
2758 printf("\n class: %d prob: %f",
2759 tree->mod->ras->parent_class_number,
2760 curr_mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number]);
2761 }
2762 }
2763
2764 if(print == YES)
2765 {
2766 PhyML_Fprintf(fp,"%4d\t%4d\t",site+1,d->num);
2767 for(i=0;i<Ns;i++)
2768 {
2769 PhyML_Fprintf(fp,"%.4f\t",p[i]);
2770 }
2771 PhyML_Fprintf(fp,"\n");
2772 fflush(NULL);
2773 }
2774
2775 tree = tree->next;
2776 d = d->next;
2777
2778
2779 }
2780 while(tree && d && tree->is_mixt_tree == NO);
2781 }
2782
2783 Free(p);
2784 curr_mixt_tree = curr_mixt_tree->next_mixt;
2785 curr_mixt_d = curr_mixt_d->next_mixt;
2786 }
2787 while(curr_mixt_tree != NULL);
2788 }
2789 }
2790
2791 //////////////////////////////////////////////////////////////
2792 //////////////////////////////////////////////////////////////
2793
2794 // First and second derivative of the log-likelihood with respect
2795 // to the length of edge b
MIXT_dLk(phydbl * l,t_edge * mixt_b,t_tree * mixt_tree)2796 phydbl MIXT_dLk(phydbl *l, t_edge *mixt_b, t_tree *mixt_tree)
2797 {
2798 t_tree *tree,*cpy_mixt_tree;
2799 t_edge *b,*cpy_mixt_b;
2800 unsigned int site, class, nclasses, ns, state;
2801 phydbl *sum_scale_left_cat,*sum_scale_rght_cat;
2802 phydbl sum,mult;
2803 phydbl *lk,*dlk;
2804 phydbl site_lk,site_dlk;
2805 phydbl log_site_lk,inv_site_lk;
2806 int num_prec_issue;
2807 phydbl r_mat_weight_sum, e_frq_weight_sum, sum_probas;
2808 phydbl len;
2809 phydbl *expl,*dot_prod;
2810 phydbl rr;
2811 phydbl ev,expevlen;
2812 phydbl one_m_pinv;
2813 short int length_found;
2814
2815 tree = NULL;
2816 b = NULL;
2817 expl = NULL;
2818 lk = NULL;
2819 dlk = NULL;
2820 cpy_mixt_tree = mixt_tree;
2821 cpy_mixt_b = mixt_b;
2822 len = -1.;
2823
2824 if(mixt_tree->update_eigen_lr == YES)
2825 {
2826 MIXT_Update_Eigen_Lr(mixt_b,mixt_tree);
2827 }
2828
2829 /*! Make sure that l is one of the lengths of mixt_b */
2830 b = mixt_b;
2831 while(b) { if(&(b->l->v) == l) break; b = b->next; }
2832 assert(b != NULL);
2833
2834
2835 do /*! Consider each element of the data partition */
2836 {
2837 b = mixt_b;
2838 tree = mixt_tree;
2839 length_found = NO;
2840 do
2841 {
2842 if(&(b->l->v) == l)
2843 {
2844 length_found = YES;
2845 break;
2846 }
2847
2848 b = b->next;
2849 tree = tree->next;
2850 }
2851 while(tree && tree->is_mixt_tree == NO);
2852
2853 /*! For computational efficiency, the dlk and lk computation for
2854 data partition elements corresponding to trees in which l is not
2855 found could be skipped. We compute dlk and lk nonetheless so that
2856 these two values are the derivative and lielihood computed for the
2857 whole data set, not juste the data partitions that "contain" l
2858 */
2859
2860 tree = mixt_tree;
2861 do
2862 {
2863 tree->c_lnL = .0;
2864 tree->c_dlnL = .0;
2865 tree = tree->next;
2866 }
2867 while(tree);
2868
2869 ns = mixt_tree->mod->ns;
2870 nclasses = MIXT_Mixt_Size(mixt_tree);
2871
2872 lk = (phydbl *)mCalloc(nclasses,sizeof(phydbl));
2873 dlk = (phydbl *)mCalloc(nclasses,sizeof(phydbl));
2874
2875 sum_scale_left_cat = (phydbl *)mCalloc(nclasses,sizeof(phydbl));
2876 sum_scale_rght_cat = (phydbl *)mCalloc(nclasses,sizeof(phydbl));
2877
2878 r_mat_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(mixt_tree->next->mod->r_mat_weight);
2879 e_frq_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(mixt_tree->next->mod->e_frq_weight);
2880 sum_probas = MIXT_Get_Sum_Of_Probas_Across_Mixtures(r_mat_weight_sum,e_frq_weight_sum,mixt_tree);
2881
2882 // Fill in expl vector for each rate class
2883 tree = mixt_tree->next;
2884 b = mixt_b->next;
2885 class = 0;
2886 do
2887 {
2888 if(tree->mod->ras->invar == YES) rr = 0.0;
2889 else
2890 {
2891 rr = 1.0;
2892 rr *= tree->mod->br_len_mult->v;
2893 rr *= mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number];
2894 }
2895
2896
2897 if(length_found == YES) len = (*l) * rr;
2898 else len = b->l->v * rr;
2899
2900 if(isinf(len) || isnan(len))
2901 {
2902 PhyML_Fprintf(stderr,"\n. len=%f rr=%f l=%f",len,rr,*l);
2903 assert(FALSE);
2904 }
2905
2906 if(tree->mod->ras->invar == NO)
2907 {
2908 if(len < tree->mod->l_min) len = tree->mod->l_min;
2909 else if(len > tree->mod->l_max) len = tree->mod->l_max;
2910 }
2911 else
2912 {
2913 len = 0.0;
2914 }
2915
2916
2917 expl = mixt_tree->expl;
2918
2919 for(state=0;state<tree->mod->ns;++state)
2920 {
2921 ev = tree->mod->eigen->e_val[state];
2922 expevlen = exp(ev*len);
2923
2924 expl[class*2*ns + 2*state] = expevlen;
2925 expl[class*2*ns + 2*state + 1] = expevlen*ev*rr;
2926 }
2927
2928 class++;
2929
2930 tree = tree->next;
2931 b = b->next;
2932 }
2933 while(tree && tree->is_mixt_tree == NO);
2934
2935 assert(class == nclasses);
2936
2937 mixt_tree->c_lnL = .0;
2938 mixt_tree->c_dlnL = .0;
2939
2940 for(site=0;site<mixt_tree->n_pattern;++site)
2941 {
2942 for(class=0;class<nclasses;++class)
2943 {
2944 dlk[class] = 0.0;
2945 lk[class] = 0.0;
2946 }
2947
2948
2949 b = mixt_b->next;
2950 tree = mixt_tree->next;
2951 class = 0;
2952 /*! For all classes in the mixture */
2953 do
2954 {
2955 if(tree->mod->ras->invar == NO && tree->data->wght[tree->curr_site] > SMALL)
2956 {
2957 tree->curr_site = site;
2958 dot_prod = tree->dot_prod + site * ns;
2959 expl = mixt_tree->expl + 2*ns*class;
2960
2961 if(tree->mod->io->datatype == NT || tree->mod->io->datatype == AA)
2962 {
2963 #if (defined(__AVX__))
2964 AVX_Lk_dLk_Core_One_Class_Eigen_Lr(dot_prod,
2965 expl ? expl : NULL,
2966 ns,lk+class,dlk+class);
2967 #elif (defined(__SSE3__))
2968 SSE_Lk_dLk_Core_One_Class_Eigen_Lr(dot_prod,
2969 expl ? expl : NULL,
2970 ns,lk+class,dlk+class);
2971 #else
2972 Lk_dLk_Core_One_Class_Eigen_Lr(dot_prod,
2973 expl ? expl : NULL,
2974 ns,lk+class,dlk+class);
2975 #endif
2976 }
2977 else
2978 {
2979 Lk_dLk_Core_One_Class_Eigen_Lr(dot_prod,
2980 expl ? expl : NULL,
2981 ns,lk+class,dlk+class);
2982 }
2983
2984 tree->apply_lk_scaling = YES;
2985
2986
2987 sum_scale_left_cat[tree->mod->ras->parent_class_number] =
2988 (b->sum_scale_left)?
2989 (b->sum_scale_left[site]):
2990 (0.0);
2991
2992 sum_scale_rght_cat[tree->mod->ras->parent_class_number] =
2993 (b->sum_scale_rght)?
2994 (b->sum_scale_rght[site]):
2995 (0.0);
2996
2997 sum =
2998 sum_scale_left_cat[tree->mod->ras->parent_class_number] +
2999 sum_scale_rght_cat[tree->mod->ras->parent_class_number];
3000
3001 if(sum > 1024.)
3002 {
3003 /* PhyML_Fprintf(stderr,"\n. Numerical precision issue detected (sum = %g)!!!",sum); */
3004 sum = 1023.;
3005 tree->mixt_tree->numerical_warning = YES;
3006 tree->numerical_warning = YES;
3007 }
3008
3009 mult = pow(2,sum);
3010
3011 /* PhyML_Printf("\n> tree: %p class: %d lk: %f dlk: %g sum: %g site: %d dot_prod[0]: %g expl[0]: %g expl[1]: %g expl[2]: %g",tree,class,log(lk[class]),dlk[class],sum,site,dot_prod[0],expl[0],expl[1],expl[2]); */
3012
3013 lk[class] /= mult;
3014 dlk[class] /= mult;
3015
3016 }
3017
3018 tree = tree->next;
3019 b = b->next;
3020 class++;
3021
3022 }
3023 while(tree && tree->is_mixt_tree == NO);
3024 // done with all trees in the mixture for this partition element.
3025
3026 tree = mixt_tree->next;
3027 class = 0;
3028 site_lk = .0;
3029 site_dlk = .0;
3030
3031 if(mixt_tree->mod->ras->invar == YES) one_m_pinv = 1. - mixt_tree->mod->ras->pinvar->v;
3032 else one_m_pinv = 1.;
3033
3034
3035 do
3036 {
3037 if(tree->mod->ras->invar == NO && mixt_tree->data->wght[tree->curr_site] > SMALL)
3038 {
3039 site_lk +=
3040 lk[class] *
3041 mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number] *
3042 tree->mod->r_mat_weight->v / r_mat_weight_sum *
3043 tree->mod->e_frq_weight->v / e_frq_weight_sum /
3044 sum_probas;
3045
3046
3047 site_dlk +=
3048 dlk[class] *
3049 one_m_pinv *
3050 mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number] *
3051 tree->mod->r_mat_weight->v / r_mat_weight_sum *
3052 tree->mod->e_frq_weight->v / e_frq_weight_sum /
3053 sum_probas;
3054
3055 /* PhyML_Printf("\n. MIXT_DLK site: %d mixt_tree: %p tree: %p class: %d l: %p lk: %g dlk: %g", */
3056 /* site, */
3057 /* mixt_tree, */
3058 /* tree, */
3059 /* class, */
3060 /* l, */
3061 /* lk[class], */
3062 /* dlk[class]); */
3063
3064 }
3065
3066 class++;
3067 tree = tree->next;
3068 }
3069 while(tree && tree->is_mixt_tree == NO);
3070
3071
3072
3073 /* Scaling for invariants */
3074 if(mixt_tree->mod->ras->invar == YES)
3075 {
3076 num_prec_issue = NO;
3077
3078 tree = mixt_tree->next;
3079 while(tree->mod->ras->invar == NO)
3080 {
3081 tree = tree->next;
3082 if(!tree || tree->is_mixt_tree == YES)
3083 {
3084 PhyML_Fprintf(stderr,"\n. tree: %p",tree);
3085 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d",__FILE__,__LINE__);
3086 Exit("\n");
3087 }
3088 }
3089
3090 tree->apply_lk_scaling = YES;
3091
3092 /*! 'tree' will give the correct state frequencies (as opposed to mixt_tree) */
3093 inv_site_lk = Invariant_Lk(0,site,&num_prec_issue,tree);
3094
3095
3096 if(num_prec_issue == YES) // inv_site_lk >> site_lk
3097 {
3098 site_lk = inv_site_lk * mixt_tree->mod->ras->pinvar->v;
3099 }
3100 else
3101 {
3102 site_lk = site_lk * (1. - mixt_tree->mod->ras->pinvar->v) + inv_site_lk * mixt_tree->mod->ras->pinvar->v;
3103 }
3104 }
3105
3106
3107 log_site_lk = log(site_lk);
3108
3109 if(isinf(log_site_lk) || isnan(log_site_lk))
3110 {
3111 PhyML_Fprintf(stderr,"\n. site = %d",site);
3112 PhyML_Fprintf(stderr,"\n. invar = %d",mixt_tree->data->invar[site]);
3113 PhyML_Fprintf(stderr,"\n. mixt = %d",mixt_tree->is_mixt_tree);
3114 PhyML_Fprintf(stderr,"\n. lk = %G log(lk) = %f < %G",site_lk,log_site_lk,-BIG);
3115 for(class=0;class<mixt_tree->mod->ras->n_catg;++class) PhyML_Printf("\n. rr=%f p=%f",mixt_tree->mod->ras->gamma_rr->v[class],mixt_tree->mod->ras->gamma_r_proba->v[class]);
3116 PhyML_Fprintf(stderr,"\n. pinv = %G",mixt_tree->mod->ras->pinvar->v);
3117 PhyML_Fprintf(stderr,"\n. bl mult = %G",mixt_tree->mod->br_len_mult->v);
3118 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d.\n",__FILE__,__LINE__);
3119 Exit("\n");
3120 }
3121
3122
3123 mixt_tree->c_lnL += mixt_tree->data->wght[site] * log_site_lk;
3124 mixt_tree->c_dlnL += mixt_tree->data->wght[site] * ( site_dlk / site_lk );
3125
3126 /* PhyML_Printf("\n. site: %4d lnL: %15f dlnL: %15f",site,site_lk,site_dlk); */
3127 }
3128
3129 Free(sum_scale_left_cat);
3130 Free(sum_scale_rght_cat);
3131
3132 Free(lk);
3133 Free(dlk);
3134
3135 mixt_tree = mixt_tree->next_mixt;
3136 mixt_b = mixt_b->next_mixt;
3137 }
3138 while(mixt_tree);
3139
3140
3141 // Sum c_dlnL across all partition element and set each c_dlnL equal to
3142 // this sum
3143
3144 sum = .0;
3145 mixt_tree = cpy_mixt_tree;
3146 do
3147 {
3148 sum += mixt_tree->c_dlnL;
3149 mixt_tree = mixt_tree->next_mixt;
3150 }
3151 while(mixt_tree);
3152
3153 mixt_tree = cpy_mixt_tree;
3154 do
3155 {
3156 mixt_tree->c_dlnL = sum;
3157 mixt_tree = mixt_tree->next_mixt;
3158 }
3159 while(mixt_tree);
3160
3161
3162
3163 // Do the same with likelihood values
3164 sum = .0;
3165 mixt_tree = cpy_mixt_tree;
3166 do
3167 {
3168 sum += mixt_tree->c_lnL;
3169 mixt_tree = mixt_tree->next_mixt;
3170 }
3171 while(mixt_tree);
3172
3173 mixt_tree = cpy_mixt_tree;
3174 do
3175 {
3176 mixt_tree->c_lnL = sum;
3177 mixt_tree = mixt_tree->next_mixt;
3178 }
3179 while(mixt_tree);
3180
3181 mixt_tree = cpy_mixt_tree;
3182 mixt_b = cpy_mixt_b;
3183
3184 return mixt_tree->c_lnL;
3185 }
3186
3187 //////////////////////////////////////////////////////////////
3188 //////////////////////////////////////////////////////////////
3189 // Returns the number of trees (including mixt_trees) in the
3190 // whole model.
MIXT_Part_Mixt_Size(t_tree * mixt_tree)3191 int MIXT_Part_Mixt_Size(t_tree *mixt_tree)
3192 {
3193 if(mixt_tree->is_mixt_tree == NO) return 1;
3194 else
3195 {
3196 int num;
3197 t_tree *tree;
3198
3199 num = 0;
3200 tree = mixt_tree;
3201 do
3202 {
3203 num++;
3204 tree = tree->next;
3205 }
3206 while(tree);
3207 return num;
3208 }
3209 return -1;
3210 }
3211
3212 //////////////////////////////////////////////////////////////
3213 //////////////////////////////////////////////////////////////
3214
3215 // Returns the number of trees in the partition element corresponding to mixt_tree.
MIXT_Mixt_Size(t_tree * mixt_tree)3216 int MIXT_Mixt_Size(t_tree *mixt_tree)
3217 {
3218 if(mixt_tree->is_mixt_tree == NO) return 1;
3219 else
3220 {
3221 int num;
3222 t_tree *tree;
3223
3224 num = 0;
3225 tree = mixt_tree->next;
3226 do
3227 {
3228 num++;
3229 tree = tree->next;
3230 }
3231 while(tree && tree->is_mixt_tree == NO);
3232
3233 return num;
3234 }
3235 return -1;
3236 }
3237
3238 //////////////////////////////////////////////////////////////
3239 //////////////////////////////////////////////////////////////
3240
MIXT_Set_Both_Sides(int yesno,t_tree * mixt_tree)3241 void MIXT_Set_Both_Sides(int yesno, t_tree *mixt_tree)
3242 {
3243 t_tree *tree;
3244
3245 assert(mixt_tree->is_mixt_tree == YES);
3246
3247 tree = mixt_tree;
3248
3249 do
3250 {
3251 if(tree->is_mixt_tree == YES) tree = tree->next;
3252 Set_Both_Sides(yesno,tree);
3253 tree = tree->next;
3254 }
3255 while(tree);
3256
3257 }
3258
3259 //////////////////////////////////////////////////////////////
3260 //////////////////////////////////////////////////////////////
3261
MIXT_Set_Model_Parameters(t_mod * mixt_mod)3262 void MIXT_Set_Model_Parameters(t_mod *mixt_mod)
3263 {
3264 t_mod *mod = mixt_mod->next;
3265
3266 if(mod != NULL)
3267 {
3268 do
3269 {
3270 Set_Model_Parameters(mod);
3271 mod = mod->next;
3272 }
3273 while(mod);
3274 }
3275 }
3276
3277 //////////////////////////////////////////////////////////////
3278 //////////////////////////////////////////////////////////////
3279
MIXT_Update_Eigen(t_mod * mixt_mod)3280 void MIXT_Update_Eigen(t_mod *mixt_mod)
3281 {
3282 t_mod *mod;
3283
3284 mod = mixt_mod;
3285
3286 do
3287 {
3288 if(mod->is_mixt_mod) mod = mod->next;
3289 if(!Update_Boundaries(mod)) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
3290 if(!Update_Eigen(mod)) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
3291 mod = mod->next;
3292 }
3293 while(mod);
3294
3295 }
3296
3297 //////////////////////////////////////////////////////////////
3298 //////////////////////////////////////////////////////////////
3299
MIXT_Update_Eigen_Lr(t_edge * mixt_b,t_tree * mixt_tree)3300 void MIXT_Update_Eigen_Lr(t_edge *mixt_b, t_tree *mixt_tree)
3301 {
3302 t_edge *b;
3303 t_tree *tree;
3304
3305 tree = mixt_tree;
3306 b = mixt_b;
3307
3308 do
3309 {
3310 if(tree->is_mixt_tree == YES)
3311 {
3312 tree = tree->next;
3313 b = b->next;
3314 }
3315
3316 Update_Eigen_Lr(b,tree);
3317
3318 tree = tree->next;
3319 b = b->next;
3320 }
3321 while(tree);
3322 }
3323
3324
3325 //////////////////////////////////////////////////////////////
3326 //////////////////////////////////////////////////////////////
3327
MIXT_Set_Update_Eigen(int yn,t_mod * mixt_mod)3328 void MIXT_Set_Update_Eigen(int yn, t_mod *mixt_mod)
3329 {
3330 t_mod *mod;
3331
3332 mod = mixt_mod;
3333 do
3334 {
3335 mod->update_eigen = yn;
3336 mod = mod->next;
3337 }
3338 while(mod);
3339 }
3340
3341 //////////////////////////////////////////////////////////////
3342 //////////////////////////////////////////////////////////////
3343
MIXT_Set_Update_Eigen_Lr(int yn,t_tree * mixt_tree)3344 void MIXT_Set_Update_Eigen_Lr(int yn, t_tree *mixt_tree)
3345 {
3346 t_tree *tree;
3347 int dum;
3348
3349 tree = mixt_tree->next;
3350 do
3351 {
3352 dum = mixt_tree->is_mixt_tree;
3353 mixt_tree->is_mixt_tree = NO;
3354
3355 Set_Update_Eigen_Lr(yn,tree);
3356
3357 mixt_tree->is_mixt_tree = dum;
3358
3359 tree = tree->next;
3360 }
3361 while(tree);
3362 }
3363
3364 //////////////////////////////////////////////////////////////
3365 //////////////////////////////////////////////////////////////
3366
MIXT_Set_Use_Eigen_Lr(int yn,t_tree * mixt_tree)3367 void MIXT_Set_Use_Eigen_Lr(int yn, t_tree *mixt_tree)
3368 {
3369 t_tree *tree;
3370 int dum;
3371
3372 tree = mixt_tree->next;
3373 do
3374 {
3375 dum = mixt_tree->is_mixt_tree;
3376 mixt_tree->is_mixt_tree = NO;
3377
3378 Set_Use_Eigen_Lr(yn,tree);
3379
3380 mixt_tree->is_mixt_tree = dum;
3381
3382 tree = tree->next;
3383 }
3384 while(tree);
3385 }
3386
3387 //////////////////////////////////////////////////////////////
3388 //////////////////////////////////////////////////////////////
3389
MIXT_Sample_Ancestral_Seq(int fullmutmap,int fromprior,t_tree * mixt_tree)3390 void MIXT_Sample_Ancestral_Seq(int fullmutmap, int fromprior, t_tree *mixt_tree)
3391 {
3392 t_tree *tree,*loc_mixt_tree;
3393 phydbl r_mat_weight_sum, e_frq_weight_sum, sum_probas, *class_proba;
3394 int i;
3395
3396 r_mat_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(mixt_tree->next->mod->r_mat_weight);
3397 e_frq_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(mixt_tree->next->mod->e_frq_weight);
3398 sum_probas = MIXT_Get_Sum_Of_Probas_Across_Mixtures(r_mat_weight_sum,e_frq_weight_sum,mixt_tree);
3399
3400
3401 tree = mixt_tree;
3402 loc_mixt_tree = tree;
3403
3404 // For each element in the partition (i.e., each gene), select a
3405 // class of the mixture according to its post prob.
3406 do
3407 {
3408 tree = loc_mixt_tree->next;
3409
3410 class_proba = NULL;
3411 i = 0;
3412 do
3413 {
3414 if(i == 0) class_proba = (phydbl *)mCalloc(1,sizeof(phydbl));
3415 else class_proba = (phydbl *)mRealloc(class_proba,i+1,sizeof(phydbl));
3416
3417 class_proba[i] =
3418 tree->site_lk_cat[0] *
3419 loc_mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number] *
3420 tree->mod->r_mat_weight->v / r_mat_weight_sum *
3421 tree->mod->e_frq_weight->v / e_frq_weight_sum /
3422 sum_probas;
3423
3424 tree = tree->next;
3425 i++;
3426 }
3427 while(tree && tree->is_mixt_tree == NO);
3428
3429 tree = loc_mixt_tree->next;
3430 i = Sample_i_With_Proba_pi(class_proba,i);
3431
3432 do
3433 {
3434 i--;
3435 if(i < 0) break;
3436 tree = tree->next;
3437 assert(tree);
3438 }
3439 while(1);
3440
3441 assert(tree->is_mixt_tree == NO);
3442
3443 Sample_Ancestral_Seq(fullmutmap,fromprior,tree);
3444
3445 Free(class_proba);
3446
3447 loc_mixt_tree = loc_mixt_tree->next_mixt;
3448 }
3449 while(loc_mixt_tree);
3450
3451 }
3452
3453 //////////////////////////////////////////////////////////////
3454 //////////////////////////////////////////////////////////////
3455 // Propagate changes in the "master" tree to other trees. Partial
3456 // likelihood vectors (and scaling stuff) is not update here so
3457 // requires full tree traversal to be updated properly.
MIXT_Propagate_Tree_Update(t_tree * mixt_tree)3458 void MIXT_Propagate_Tree_Update(t_tree *mixt_tree)
3459 {
3460 t_tree *tree;
3461 int i;
3462
3463 assert(!mixt_tree->prev);
3464
3465 tree = mixt_tree->next;
3466 if(!tree) return;
3467
3468 do
3469 {
3470 for(i=0;i<2*mixt_tree->n_otu-1;++i)
3471 {
3472 tree->a_nodes[i]->v[0] = mixt_tree->a_nodes[i]->v[0] ? tree->a_nodes[mixt_tree->a_nodes[i]->v[0]->num] : NULL;
3473 tree->a_nodes[i]->v[1] = mixt_tree->a_nodes[i]->v[1] ? tree->a_nodes[mixt_tree->a_nodes[i]->v[1]->num] : NULL;
3474 tree->a_nodes[i]->v[2] = mixt_tree->a_nodes[i]->v[2] ? tree->a_nodes[mixt_tree->a_nodes[i]->v[2]->num] : NULL;
3475
3476 if(tree->times != NULL) tree->times->nd_t[i] = mixt_tree->times->nd_t[i];
3477 }
3478
3479 Connect_Edges_To_Nodes_Serial(tree);
3480
3481 Reorganize_Edges_Given_Lk_Struct(tree);
3482 Init_Partial_Lk_Tips_Double(tree);
3483
3484 if(mixt_tree->n_root != NULL)
3485 {
3486 assert(mixt_tree->e_root);
3487 assert(mixt_tree->n_root->v[1]);
3488 assert(mixt_tree->n_root->v[2]);
3489 assert(mixt_tree->n_root->b[1]);
3490 assert(mixt_tree->n_root->b[2]);
3491
3492 tree->n_root = tree->a_nodes[mixt_tree->n_root->num];
3493 tree->e_root = tree->a_edges[mixt_tree->e_root->num];
3494
3495 tree->n_root->v[1] = tree->a_nodes[mixt_tree->n_root->v[1]->num];
3496 tree->n_root->v[2] = tree->a_nodes[mixt_tree->n_root->v[2]->num];
3497
3498 tree->n_root->b[1] = tree->a_edges[mixt_tree->n_root->b[1]->num];
3499 tree->n_root->b[2] = tree->a_edges[mixt_tree->n_root->b[2]->num];
3500
3501 tree->n_root->b[1]->left = tree->n_root;
3502 tree->n_root->b[2]->left = tree->n_root;
3503 tree->n_root->b[1]->rght = tree->n_root->v[1];
3504 tree->n_root->b[2]->rght = tree->n_root->v[2];
3505 }
3506
3507 Update_Ancestors(tree->n_root,tree->n_root->v[2],tree);
3508 Update_Ancestors(tree->n_root,tree->n_root->v[1],tree);
3509
3510 tree = tree->next;
3511
3512 }
3513 while(tree);
3514 }
3515
3516 //////////////////////////////////////////////////////////////
3517 //////////////////////////////////////////////////////////////
3518
MIXT_Set_Bl_From_Rt(int yn,t_tree * mixt_tree)3519 void MIXT_Set_Bl_From_Rt(int yn, t_tree *mixt_tree)
3520 {
3521 t_tree *tree;
3522
3523 tree = mixt_tree;
3524 do
3525 {
3526 assert(tree->rates);
3527 tree->rates->bl_from_rt = yn;
3528 tree = tree->next;
3529 }
3530 while(tree);
3531 }
3532
3533 //////////////////////////////////////////////////////////////
3534 //////////////////////////////////////////////////////////////
3535
MIXT_Copy_Tree(t_tree * ori,t_tree * cpy)3536 void MIXT_Copy_Tree(t_tree *ori, t_tree *cpy)
3537 {
3538 int ori_is_mixt_tree,cpy_is_mixt_tree;
3539
3540 assert(!((cpy && !ori) || (!cpy && ori)));
3541
3542 if(cpy == NULL || ori == NULL) return;
3543
3544 if(ori->is_mixt_tree == YES && cpy->is_mixt_tree == YES)
3545 {
3546 do
3547 {
3548 ori_is_mixt_tree = ori->is_mixt_tree;
3549 ori->is_mixt_tree = NO;
3550 cpy_is_mixt_tree = cpy->is_mixt_tree;
3551 cpy->is_mixt_tree = NO;
3552
3553 Copy_Tree(ori,cpy);
3554
3555 cpy->is_mixt_tree = cpy_is_mixt_tree;
3556 ori->is_mixt_tree = ori_is_mixt_tree;
3557
3558 ori = ori->next;
3559 cpy = cpy->next;
3560 }
3561 while(cpy);
3562 }
3563 else if(ori->is_mixt_tree == NO && cpy->is_mixt_tree == YES)
3564 {
3565 do
3566 {
3567 cpy_is_mixt_tree = cpy->is_mixt_tree;
3568 cpy->is_mixt_tree = NO;
3569
3570 Copy_Tree(ori,cpy);
3571
3572 cpy->is_mixt_tree = cpy_is_mixt_tree;
3573
3574 cpy = cpy->next;
3575 }
3576 while(cpy);
3577 }
3578 else if(ori->is_mixt_tree == YES && cpy->is_mixt_tree == NO)
3579 {
3580 ori_is_mixt_tree = ori->is_mixt_tree;
3581 ori->is_mixt_tree = NO;
3582 Copy_Tree(ori,cpy);
3583 ori->is_mixt_tree = ori_is_mixt_tree;
3584 }
3585 }
3586
3587 //////////////////////////////////////////////////////////////
3588 //////////////////////////////////////////////////////////////
3589
MIXT_Init_NNI_Score(phydbl val,t_edge * mixt_b,t_tree * mixt_tree)3590 void MIXT_Init_NNI_Score(phydbl val, t_edge *mixt_b, t_tree *mixt_tree)
3591 {
3592 t_edge *b;
3593 t_tree *tree;
3594
3595 tree = mixt_tree->next;
3596 b = mixt_b->next;
3597 do
3598 {
3599 if(tree->is_mixt_tree)
3600 {
3601 tree = tree->next;
3602 b = b->next;
3603 }
3604
3605 Init_NNI_Score(val,b,tree);
3606 tree = tree->next;
3607 b = b->next;
3608 }
3609 while(tree);
3610
3611 }
3612
3613 //////////////////////////////////////////////////////////////
3614 //////////////////////////////////////////////////////////////
3615
MIXT_Duplicate_Tree(t_tree * ori)3616 t_tree *MIXT_Duplicate_Tree(t_tree *ori)
3617 {
3618 t_tree *cpy,*first,*prev;
3619 int dum;
3620
3621 ori->is_mixt_tree = NO;
3622 first = Duplicate_Tree(ori);
3623 first->prev = NULL;
3624 first->is_mixt_tree = YES;
3625 ori->is_mixt_tree = YES;
3626
3627 ori = ori->next;
3628 prev = cpy = first;
3629 do
3630 {
3631 dum = ori->is_mixt_tree;
3632 ori->is_mixt_tree = NO;
3633
3634 prev = cpy;
3635
3636 cpy = Duplicate_Tree(ori);
3637
3638 ori->is_mixt_tree = dum;
3639 cpy->is_mixt_tree = ori->is_mixt_tree;
3640
3641 cpy->prev = prev;
3642 prev->next = cpy;
3643
3644 ori = ori->next;
3645 }
3646 while(ori);
3647
3648 cpy->next = NULL;
3649
3650 return(first);
3651 }
3652
3653 //////////////////////////////////////////////////////////////
3654 //////////////////////////////////////////////////////////////
3655
MIXT_Print_Site_Lk(t_tree * mixt_tree,FILE * fp)3656 void MIXT_Print_Site_Lk(t_tree *mixt_tree, FILE *fp)
3657 {
3658 t_tree *tree;
3659 int site,catg;
3660 char *s;
3661 phydbl postmean,sum;
3662
3663 assert(mixt_tree->is_mixt_tree == YES);
3664 assert(mixt_tree->io->print_site_lnl == YES);
3665
3666
3667 if(!mixt_tree->io->print_trace)
3668 {
3669 tree = mixt_tree;
3670 do
3671 {
3672 PhyML_Fprintf(fp,"Note : P(D|M) is the probability of site D given the model M (i.e., the site likelihood)\n");
3673 if(tree->mod->ras->n_catg > 1 || tree->mod->ras->invar)
3674 {
3675 PhyML_Fprintf(fp,"P*(D|M,rr[x]) is the scaled probability of site D given the model M and the relative rate\n");
3676 PhyML_Fprintf(fp,"of evolution rr[x], where x is the class of rate to be considered.\n");
3677 PhyML_Fprintf(fp,"The actual conditional probability is given by P*(D|M,rr[x])/2^F, where\n");
3678 PhyML_Fprintf(fp,"F is the scaling factor (see column 'Scaler').\n");
3679 PhyML_Fprintf(fp,"For invariant sites, P(D|M,rr[0]=0) is the actual conditional probability\n");
3680 PhyML_Fprintf(fp,"(i.e., it is not scaled).\n");
3681 break;
3682 }
3683 tree = tree->next_mixt;
3684 }
3685 while(tree);
3686 PhyML_Fprintf(fp,"\n");
3687
3688
3689 s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3690
3691 do
3692 {
3693 PhyML_Fprintf(fp,"Alignment file name: %s\n\n",mixt_tree->io->in_align_file);
3694
3695 sprintf(s,"Site");
3696 PhyML_Fprintf(fp, "%-12s",s);
3697
3698 sprintf(s,"P(D|M)");
3699 PhyML_Fprintf(fp,"%-15s",s);
3700
3701 sprintf(s,"Scaler");
3702 PhyML_Fprintf(fp,"%-7s",s);
3703
3704 sprintf(s,"Pattern");
3705 PhyML_Fprintf(fp, "%-9s",s);
3706
3707 if(mixt_tree->mod->ras->n_catg > 1)
3708 {
3709 for(catg=0;catg<mixt_tree->mod->ras->n_catg;catg++)
3710 {
3711 sprintf(s,"P*(D|M,rr[%d]=%5.4f)",catg+1,mixt_tree->mod->ras->gamma_rr->v[catg]);
3712 PhyML_Fprintf(fp,"%-23s",s);
3713 }
3714
3715 sprintf(s,"Posterior mean");
3716 PhyML_Fprintf(fp,"%-22s",s);
3717 }
3718
3719
3720 if(mixt_tree->mod->ras->invar)
3721 {
3722 sprintf(s,"P(D|M,rr[0]=0)");
3723 PhyML_Fprintf(fp,"%-16s",s);
3724 }
3725
3726 sprintf(s,"NDistinctStates");
3727 PhyML_Fprintf(fp,"%-16s",s);
3728
3729 PhyML_Fprintf(fp,"\n");
3730
3731 assert(mixt_tree->next->is_mixt_tree == NO);
3732 Init_Ui_Tips(mixt_tree->next);
3733
3734 for(site=0;site<mixt_tree->data->init_len;site++)
3735 {
3736 PhyML_Fprintf(fp,"%-12d",site+1);
3737
3738 PhyML_Fprintf(fp,"%-15g",mixt_tree->cur_site_lk[mixt_tree->data->sitepatt[site]]);
3739
3740 PhyML_Fprintf(fp,"%-7d",mixt_tree->fact_sum_scale[mixt_tree->data->sitepatt[site]]);
3741
3742 PhyML_Fprintf(fp,"%-9d",mixt_tree->data->sitepatt[site]);
3743
3744 if(mixt_tree->mod->ras->n_catg > 1)
3745 {
3746 tree = mixt_tree->next;
3747 do
3748 {
3749 PhyML_Fprintf(fp,"%-23g",tree->unscaled_site_lk_cat[tree->data->sitepatt[site]]);
3750 tree = tree->next;
3751 }
3752 while(tree && tree->is_mixt_tree == NO);
3753
3754
3755 tree = mixt_tree->next;
3756 postmean = .0;
3757 sum = .0;
3758 do
3759 {
3760 postmean +=
3761 mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number] *
3762 tree->unscaled_site_lk_cat[tree->data->sitepatt[site]] *
3763 mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number];
3764
3765 sum +=
3766 tree->unscaled_site_lk_cat[tree->data->sitepatt[site]] *
3767 mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number];
3768
3769 tree = tree->next;
3770 }
3771 while(tree && tree->is_mixt_tree == NO);
3772
3773 postmean /= sum;
3774
3775 PhyML_Fprintf(fp,"%-22g",postmean);
3776 }
3777
3778 if(mixt_tree->mod->ras->invar)
3779 {
3780 if((phydbl)mixt_tree->data->invar[mixt_tree->data->sitepatt[site]] > -0.5)
3781 PhyML_Fprintf(fp,"%-16g",mixt_tree->mod->e_frq->pi->v[mixt_tree->data->invar[mixt_tree->data->sitepatt[site]]]);
3782 else
3783 PhyML_Fprintf(fp,"%-16g",0.0);
3784 }
3785
3786 assert(mixt_tree->next != NULL);
3787 PhyML_Fprintf(fp,"%-16d",Number_Of_Diff_States_One_Site(mixt_tree->data->sitepatt[site],mixt_tree->next));
3788
3789 PhyML_Fprintf(fp,"\n");
3790 }
3791 Free(s);
3792
3793 mixt_tree = mixt_tree->next_mixt;
3794 }
3795 while(mixt_tree);
3796 }
3797 }
3798
3799 //////////////////////////////////////////////////////////////
3800 //////////////////////////////////////////////////////////////
3801 //////////////////////////////////////////////////////////////
3802 //////////////////////////////////////////////////////////////
3803 //////////////////////////////////////////////////////////////
3804 //////////////////////////////////////////////////////////////
3805
3806
3807
3808