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 "assert.h"
14 #include "lk.h"
15 #ifdef BEAGLE
16 #include "beagle_utils.h"
17 #endif
18 #include <stdint.h>
19
20
21
22
23 //////////////////////////////////////////////////////////////
24 //////////////////////////////////////////////////////////////
25
Init_Tips_At_One_Site_Nucleotides_Float(char state,int pos,phydbl * p_lk)26 void Init_Tips_At_One_Site_Nucleotides_Float(char state, int pos, phydbl *p_lk)
27 {
28 switch(state)
29 {
30 case 'A' : p_lk[pos+0]=1.; p_lk[pos+1]=p_lk[pos+2]=p_lk[pos+3]=.0;
31 break;
32 case 'C' : p_lk[pos+1]=1.; p_lk[pos+0]=p_lk[pos+2]=p_lk[pos+3]=.0;
33 break;
34 case 'G' : p_lk[pos+2]=1.; p_lk[pos+1]=p_lk[pos+0]=p_lk[pos+3]=.0;
35 break;
36 case 'T' : p_lk[pos+3]=1.; p_lk[pos+1]=p_lk[pos+2]=p_lk[pos+0]=.0;
37 break;
38 case 'U' : p_lk[pos+3]=1.; p_lk[pos+1]=p_lk[pos+2]=p_lk[pos+0]=.0;
39 break;
40 case 'M' : p_lk[pos+0]=p_lk[pos+1]=1.; p_lk[pos+2]=p_lk[pos+3]=.0;
41 break;
42 case 'R' : p_lk[pos+0]=p_lk[pos+2]=1.; p_lk[pos+1]=p_lk[pos+3]=.0;
43 break;
44 case 'W' : p_lk[pos+0]=p_lk[pos+3]=1.; p_lk[pos+1]=p_lk[pos+2]=.0;
45 break;
46 case 'S' : p_lk[pos+1]=p_lk[pos+2]=1.; p_lk[pos+0]=p_lk[pos+3]=.0;
47 break;
48 case 'Y' : p_lk[pos+1]=p_lk[pos+3]=1.; p_lk[pos+0]=p_lk[pos+2]=.0;
49 break;
50 case 'K' : p_lk[pos+2]=p_lk[pos+3]=1.; p_lk[pos+0]=p_lk[pos+1]=.0;
51 break;
52 case 'B' : p_lk[pos+1]=p_lk[pos+2]=p_lk[pos+3]=1.; p_lk[pos+0]=.0;
53 break;
54 case 'D' : p_lk[pos+0]=p_lk[pos+2]=p_lk[pos+3]=1.; p_lk[pos+1]=.0;
55 break;
56 case 'H' : p_lk[pos+0]=p_lk[pos+1]=p_lk[pos+3]=1.; p_lk[pos+2]=.0;
57 break;
58 case 'V' : p_lk[pos+0]=p_lk[pos+1]=p_lk[pos+2]=1.; p_lk[pos+3]=.0;
59 break;
60 case 'N' : case 'X' : case '?' : case 'O' : case '-' :
61 p_lk[pos+0]=p_lk[pos+1]=p_lk[pos+2]=p_lk[pos+3]=1.;break;
62 default :
63 {
64 PhyML_Fprintf(stderr,"\n. Unknown character state : '%c'.\n",state);
65 Exit("\n. Init failed (data type supposed to be DNA)\n");
66 break;
67 }
68 }
69 }
70
71 //////////////////////////////////////////////////////////////
72 //////////////////////////////////////////////////////////////
73
Init_Tips_At_One_Site_Nucleotides_Int(char state,int pos,short int * p_pars)74 void Init_Tips_At_One_Site_Nucleotides_Int(char state, int pos, short int *p_pars)
75 {
76 switch(state)
77 {
78 case 'A' : p_pars[pos+0]=1; p_pars[pos+1]=p_pars[pos+2]=p_pars[pos+3]=0;
79 break;
80 case 'C' : p_pars[pos+1]=1; p_pars[pos+0]=p_pars[pos+2]=p_pars[pos+3]=0;
81 break;
82 case 'G' : p_pars[pos+2]=1; p_pars[pos+1]=p_pars[pos+0]=p_pars[pos+3]=0;
83 break;
84 case 'T' : p_pars[pos+3]=1; p_pars[pos+1]=p_pars[pos+2]=p_pars[pos+0]=0;
85 break;
86 case 'U' : p_pars[pos+3]=1; p_pars[pos+1]=p_pars[pos+2]=p_pars[pos+0]=0;
87 break;
88 case 'M' : p_pars[pos+0]=p_pars[pos+1]=1; p_pars[pos+2]=p_pars[pos+3]=0;
89 break;
90 case 'R' : p_pars[pos+0]=p_pars[pos+2]=1; p_pars[pos+1]=p_pars[pos+3]=0;
91 break;
92 case 'W' : p_pars[pos+0]=p_pars[pos+3]=1; p_pars[pos+1]=p_pars[pos+2]=0;
93 break;
94 case 'S' : p_pars[pos+1]=p_pars[pos+2]=1; p_pars[pos+0]=p_pars[pos+3]=0;
95 break;
96 case 'Y' : p_pars[pos+1]=p_pars[pos+3]=1; p_pars[pos+0]=p_pars[pos+2]=0;
97 break;
98 case 'K' : p_pars[pos+2]=p_pars[pos+3]=1; p_pars[pos+0]=p_pars[pos+1]=0;
99 break;
100 case 'B' : p_pars[pos+1]=p_pars[pos+2]=p_pars[pos+3]=1; p_pars[pos+0]=0;
101 break;
102 case 'D' : p_pars[pos+0]=p_pars[pos+2]=p_pars[pos+3]=1; p_pars[pos+1]=0;
103 break;
104 case 'H' : p_pars[pos+0]=p_pars[pos+1]=p_pars[pos+3]=1; p_pars[pos+2]=0;
105 break;
106 case 'V' : p_pars[pos+0]=p_pars[pos+1]=p_pars[pos+2]=1; p_pars[pos+3]=0;
107 break;
108 case 'N' : case 'X' : case '?' : case 'O' : case '-' :
109 p_pars[pos+0]=p_pars[pos+1]=p_pars[pos+2]=p_pars[pos+3]=1;break;
110 default :
111 {
112 PhyML_Fprintf(stderr,"\n. Unknown character state : '%c'.\n",state);
113 Exit("\n. Init failed (data type supposed to be DNA)\n");
114 break;
115 }
116 }
117 }
118
119 //////////////////////////////////////////////////////////////
120 //////////////////////////////////////////////////////////////
121
Init_Tips_At_One_Site_AA_Float(char aa,int pos,phydbl * p_lk)122 void Init_Tips_At_One_Site_AA_Float(char aa, int pos, phydbl *p_lk)
123 {
124 int i;
125
126 for(i=0;i<20;i++) p_lk[pos+i] = .0;
127
128 switch(aa){
129 case 'A' : p_lk[pos+0]= 1.; break;/* Alanine */
130 case 'R' : p_lk[pos+1]= 1.; break;/* Arginine */
131 case 'N' : p_lk[pos+2]= 1.; break;/* Asparagine */
132 case 'D' : p_lk[pos+3]= 1.; break;/* Aspartic acid */
133 case 'C' : p_lk[pos+4]= 1.; break;/* Cysteine */
134 case 'Q' : p_lk[pos+5]= 1.; break;/* Glutamine */
135 case 'E' : p_lk[pos+6]= 1.; break;/* Glutamic acid */
136 case 'G' : p_lk[pos+7]= 1.; break;/* Glycine */
137 case 'H' : p_lk[pos+8]= 1.; break;/* Histidine */
138 case 'I' : p_lk[pos+9]= 1.; break;/* Isoleucine */
139 case 'L' : p_lk[pos+10]=1.; break;/* Leucine */
140 case 'K' : p_lk[pos+11]=1.; break;/* Lysine */
141 case 'M' : p_lk[pos+12]=1.; break;/* Methionine */
142 case 'F' : p_lk[pos+13]=1.; break;/* Phenylalanin */
143 case 'P' : p_lk[pos+14]=1.; break;/* Proline */
144 case 'S' : p_lk[pos+15]=1.; break;/* Serine */
145 case 'T' : p_lk[pos+16]=1.; break;/* Threonine */
146 case 'W' : p_lk[pos+17]=1.; break;/* Tryptophan */
147 case 'Y' : p_lk[pos+18]=1.; break;/* Tyrosine */
148 case 'V' : p_lk[pos+19]=1.; break;/* Valine */
149
150 case 'B' : p_lk[pos+2]= 1.; break;/* Asparagine */
151 case 'Z' : p_lk[pos+5]= 1.; break;/* Glutamine */
152
153 case 'X' : case '?' : case '-' : for(i=0;i<20;i++) p_lk[pos+i] = 1.; break;
154 default :
155 {
156 PhyML_Fprintf(stderr,"\n. Unknown character state : '%c'.\n",aa);
157 Exit("\n. Init failed (data type supposed to be amino-acids)\n");
158 break;
159 }
160 }
161 }
162
163 //////////////////////////////////////////////////////////////
164 //////////////////////////////////////////////////////////////
165
Init_Tips_At_One_Site_AA_Int(char aa,int pos,short int * p_pars)166 void Init_Tips_At_One_Site_AA_Int(char aa, int pos, short int *p_pars)
167 {
168 int i;
169
170 for(i=0;i<20;i++) p_pars[pos+i] = .0;
171
172 switch(aa){
173 case 'A' : p_pars[pos+0] = 1; break;/* Alanine */
174 case 'R' : p_pars[pos+1] = 1; break;/* Arginine */
175 case 'N' : p_pars[pos+2] = 1; break;/* Asparagine */
176 case 'D' : p_pars[pos+3] = 1; break;/* Aspartic acid */
177 case 'C' : p_pars[pos+4] = 1; break;/* Cysteine */
178 case 'Q' : p_pars[pos+5] = 1; break;/* Glutamine */
179 case 'E' : p_pars[pos+6] = 1; break;/* Glutamic acid */
180 case 'G' : p_pars[pos+7] = 1; break;/* Glycine */
181 case 'H' : p_pars[pos+8] = 1; break;/* Histidine */
182 case 'I' : p_pars[pos+9] = 1; break;/* Isoleucine */
183 case 'L' : p_pars[pos+10] = 1; break;/* Leucine */
184 case 'K' : p_pars[pos+11] = 1; break;/* Lysine */
185 case 'M' : p_pars[pos+12] = 1; break;/* Methionine */
186 case 'F' : p_pars[pos+13] = 1; break;/* Phenylalanin */
187 case 'P' : p_pars[pos+14] = 1; break;/* Proline */
188 case 'S' : p_pars[pos+15] = 1; break;/* Serine */
189 case 'T' : p_pars[pos+16] = 1; break;/* Threonine */
190 case 'W' : p_pars[pos+17] = 1; break;/* Tryptophan */
191 case 'Y' : p_pars[pos+18] = 1; break;/* Tyrosine */
192 case 'V' : p_pars[pos+19] = 1; break;/* Valine */
193
194 case 'B' : p_pars[pos+2] = 1; break;/* Asparagine */
195 case 'Z' : p_pars[pos+5] = 1; break;/* Glutamine */
196
197 case 'X' : case '?' : case '-' : for(i=0;i<20;i++) p_pars[pos+i] = 1; break;
198 default :
199 {
200 PhyML_Fprintf(stderr,"\n. Unknown character state : '%c'.\n",aa);
201 Exit("\n. Init failed (data type supposed to be amino-acids)\n");
202 break;
203 }
204 }
205 }
206
207 //////////////////////////////////////////////////////////////
208 //////////////////////////////////////////////////////////////
209
Init_Tips_At_One_Site_Generic_Float(char * state,int ns,int state_len,int pos,phydbl * p_lk)210 void Init_Tips_At_One_Site_Generic_Float(char *state, int ns, int state_len, int pos, phydbl *p_lk)
211 {
212 int i;
213 int state_int;
214
215 for(i=0;i<ns;i++) p_lk[pos+i] = 0.;
216
217 if(Is_Ambigu(state,GENERIC,state_len)) for(i=0;i<ns;i++) p_lk[pos+i] = 1.;
218 else
219 {
220 char format[6];
221 sprintf(format,"%%%dd",state_len);
222 if(!sscanf(state,format,&state_int))
223 {
224 PhyML_Fprintf(stderr,"\n. state='%c'",state);
225 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d (function '%s')\n",__FILE__,__LINE__);
226 Warn_And_Exit("");
227 }
228 if(state_int > ns)
229 {
230 PhyML_Fprintf(stderr,"\n. %s %d cstate: %.2s istate: %d state_len: %d.\n",__FILE__,__LINE__,state,state_int,state_len);
231 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
232 Warn_And_Exit("");
233 }
234 p_lk[pos+state_int] = 1.;
235 /* PhyML_Printf("\n. %s %d cstate: %.2s istate: %d state_len: %d ns: %d pos: %d",__FILE__,__LINE__,state,state_int,state_len,ns,pos); */
236 }
237 }
238
239 //////////////////////////////////////////////////////////////
240 //////////////////////////////////////////////////////////////
241
Init_Tips_At_One_Site_Generic_Int(char * state,int ns,int state_len,int pos,short int * p_pars)242 void Init_Tips_At_One_Site_Generic_Int(char *state, int ns, int state_len, int pos, short int *p_pars)
243 {
244 int i;
245 int state_int;
246
247 for(i=0;i<ns;i++) p_pars[pos+i] = 0;
248
249 if(Is_Ambigu(state,GENERIC,state_len)) for(i=0;i<ns;i++) p_pars[pos+i] = 1;
250 else
251 {
252 char format[6];
253 sprintf(format,"%%%dd",state_len);
254 if(!sscanf(state,format,&state_int))
255 {
256 PhyML_Fprintf(stderr,"\n. state='%c'",state);
257 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
258 Warn_And_Exit("");
259 }
260 if(state_int > ns)
261 {
262 PhyML_Fprintf(stderr,"\n. %s %d cstate: %.2s istate: %d state_len: %d.\n",__FILE__,__LINE__,state,state_int,state_len);
263 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
264 Warn_And_Exit("");
265 }
266 p_pars[pos+state_int] = 1;
267 /* PhyML_Printf("\n* %s %d cstate: %.2s istate: %d state_len: %d ns: %d pos: %d",__FILE__,__LINE__,state,state_int,state_len,ns,pos); */
268 }
269 }
270
271 //////////////////////////////////////////////////////////////
272 //////////////////////////////////////////////////////////////
273
Get_All_Partial_Lk_Scale(t_tree * tree,t_edge * b_fcus,t_node * a,t_node * d)274 void Get_All_Partial_Lk_Scale(t_tree *tree, t_edge *b_fcus, t_node *a, t_node *d)
275 {
276 Update_Partial_Lk(tree,b_fcus,d);
277 }
278
279 //////////////////////////////////////////////////////////////
280 //////////////////////////////////////////////////////////////
281
Post_Order_Lk(t_node * a,t_node * d,t_tree * tree)282 void Post_Order_Lk(t_node *a, t_node *d, t_tree *tree)
283 {
284 int i,dir;
285
286 dir = -1;
287
288 /* PhyML_Printf("\n. %p a: %3d d: %3d root: %3d root->v1: %3d root->v2: %3d [%3d %3d %3d] (%3d %3d) (%3d %3d)", */
289 /* tree, */
290 /* a->num, */
291 /* d->num, */
292 /* tree->n_root?tree->n_root->num:-1, */
293 /* tree->n_root->v[1]->num, */
294 /* tree->n_root->v[2]->num, */
295 /* d->v[0]?d->v[0]->num:-1, */
296 /* d->v[1]?d->v[1]->num:-1, */
297 /* d->v[2]?d->v[2]->num:-1, */
298 /* tree->a_nodes[5]->num, */
299 /* tree->a_nodes[5]->v[0]->num, */
300 /* tree->mixt_tree? tree->mixt_tree->a_nodes[5]->num : -1, */
301 /* tree->mixt_tree? tree->mixt_tree->a_nodes[5]->v[0]->num : -1); */
302
303 if(d->tax) return;
304 else
305 {
306 if(tree->is_mixt_tree)
307 {
308 MIXT_Post_Order_Lk(a,d,tree);
309 return;
310 }
311
312 if(tree->n_root)
313 {
314 for(i=0;i<3;i++)
315 {
316 if(d->v[i] != a && d->b[i] != tree->e_root)
317 Post_Order_Lk(d,d->v[i],tree);
318 else dir = i;
319 }
320 }
321 else
322 {
323 for(i=0;i<3;i++)
324 {
325 if(d->v[i] != a)
326 Post_Order_Lk(d,d->v[i],tree);
327 else dir = i;
328 }
329 }
330
331 if(dir < 0)
332 {
333 PhyML_Printf("\n. a->num: %d d->num: %d d->v[0]->num: %d d->v[1]->num: %d d->v[2]->num: %d d->b[0]->num: %d d->b[1]->num: %d d->b[2]->num: %d root ? %d e_root ? %d\n",
334 a?a->num:-1,
335 d?d->num:1,
336 d->v[0]?d->v[0]->num:-1,
337 d->v[1]?d->v[1]->num:-1,
338 d->v[2]?d->v[2]->num:-1,
339 d->b[0]?d->b[0]->num:-1,
340 d->b[1]?d->b[1]->num:-1,
341 d->b[2]?d->b[2]->num:-1,
342 tree->n_root?tree->n_root->num:-1,
343 tree->e_root?tree->e_root->num:-1);
344 assert(FALSE);
345 }
346
347 if(tree->ignore_root == NO && d->b[dir] == tree->e_root)
348 {
349 if(d == tree->n_root->v[1]) Update_Partial_Lk(tree,tree->n_root->b[1],d);
350 else Update_Partial_Lk(tree,tree->n_root->b[2],d);
351 }
352 else
353 {
354 Update_Partial_Lk(tree,d->b[dir],d);
355 }
356 }
357 }
358
359 //////////////////////////////////////////////////////////////
360 //////////////////////////////////////////////////////////////
361
Pre_Order_Lk(t_node * a,t_node * d,t_tree * tree)362 void Pre_Order_Lk(t_node *a, t_node *d, t_tree *tree)
363 {
364 int i;
365
366 if(d->tax) return;
367 else
368 {
369 if(tree->is_mixt_tree)
370 {
371 MIXT_Pre_Order_Lk(a,d,tree);
372 return;
373 }
374
375 if(tree->n_root)
376 {
377 for(i=0;i<3;++i)
378 {
379 if(d->v[i] != a && d->b[i] != tree->e_root)
380 {
381 Update_Partial_Lk(tree,d->b[i],d);
382 Pre_Order_Lk(d,d->v[i],tree);
383 }
384 }
385 }
386 else
387 {
388 for(i=0;i<3;++i)
389 {
390 if(d->v[i] != a)
391 {
392 Update_Partial_Lk(tree,d->b[i],d);
393 Pre_Order_Lk(d,d->v[i],tree);
394 }
395 }
396 }
397 }
398 }
399
400 //////////////////////////////////////////////////////////////
401 //////////////////////////////////////////////////////////////
402
403 // Updates all partial likelihood vectors. Depending on whether
404 // both_sides = YES or NO, only 'up' or 'up'&'down' partials will
405 // be updates
Update_All_Partial_Lk(t_tree * tree)406 void Update_All_Partial_Lk(t_tree *tree)
407 {
408 if(tree->n_root)
409 {
410 if(tree->ignore_root == NO)
411 {
412 Post_Order_Lk(tree->n_root,tree->n_root->v[1],tree);
413 Post_Order_Lk(tree->n_root,tree->n_root->v[2],tree);
414
415 Update_Partial_Lk(tree,tree->n_root->b[1],tree->n_root);
416 Update_Partial_Lk(tree,tree->n_root->b[2],tree->n_root);
417
418 if(tree->both_sides == YES)
419 {
420 Pre_Order_Lk(tree->n_root,tree->n_root->v[2],tree);
421 Pre_Order_Lk(tree->n_root,tree->n_root->v[1],tree);
422 }
423 }
424 else
425 {
426
427 Post_Order_Lk(tree->e_root->rght,tree->e_root->left,tree);
428 Post_Order_Lk(tree->e_root->left,tree->e_root->rght,tree);
429
430 if(tree->both_sides == YES)
431 {
432 Pre_Order_Lk(tree->e_root->rght,tree->e_root->left,tree);
433 Pre_Order_Lk(tree->e_root->left,tree->e_root->rght,tree);
434 }
435 }
436 }
437 else
438 {
439 Post_Order_Lk(tree->a_nodes[tree->tip_root],tree->a_nodes[tree->tip_root]->v[0],tree);
440 if(tree->both_sides == YES)
441 Pre_Order_Lk(tree->a_nodes[tree->tip_root],tree->a_nodes[tree->tip_root]->v[0],tree);
442 }
443 }
444
445 //////////////////////////////////////////////////////////////
446 //////////////////////////////////////////////////////////////
447
Lk(t_edge * b,t_tree * tree)448 phydbl Lk(t_edge *b, t_tree *tree)
449 {
450 unsigned int br,catg,state,ambiguity_check,site;
451 phydbl len,*expl,*dot_prod,*p_lk_left,*p_lk_rght;
452
453 const unsigned int ns = tree->mod->ns;
454 const unsigned int ncatg = tree->mod->ras->n_catg;
455 const unsigned int npatterns = tree->n_pattern;
456 const unsigned int nsncatg = ns * ncatg;
457
458 tree->numerical_warning = NO;
459
460 /* if(tree->eval_alnL == NO) return UNLIKELY; */
461
462 if(b == NULL && tree->mod->s_opt->curr_opt_free_rates == YES)
463 {
464 tree->mod->s_opt->curr_opt_free_rates = NO;
465 Optimize_Free_Rate_Weights(tree,YES,YES);
466 tree->mod->s_opt->curr_opt_free_rates = YES;
467 }
468
469 #ifdef PHYREX
470 PHYREX_Ldsk_To_Tree(tree);
471 #endif
472
473
474 if(tree->is_mixt_tree == YES)
475 {
476 #ifdef BEAGLE
477 Warn_And_Exit(TODO_BEAGLE);
478 #endif
479 return MIXT_Lk(b,tree);
480 }
481
482 tree->old_lnL = tree->c_lnL;
483
484
485
486 #if (defined PHYTIME || defined INVITEE || defined PHYREX)
487 if((tree->rates) && (tree->rates->bl_from_rt)) RATES_Update_Cur_Bl(tree);
488 #endif
489
490 if(tree->rates && tree->io && tree->io->lk_approx == NORMAL)
491 {
492 #ifdef BEAGLE
493 Warn_And_Exit(TODO_BEAGLE);
494 #endif
495 tree->c_lnL = Lk_Normal_Approx(tree);
496 return tree->c_lnL;
497 }
498
499 expl = tree->expl;
500 dot_prod = tree->dot_prod;
501
502
503 if(b == NULL)
504 {
505 Update_Boundaries(tree->mod);
506 Update_RAS(tree->mod);
507 Update_Efrq(tree->mod);
508 Update_Eigen(tree->mod);
509 }
510
511
512
513 if(tree->mod->s_opt->skip_tree_traversal == NO)
514 {
515 if(!b) //Update PMat for all edges
516 {
517 for(br=0;br<2*tree->n_otu-3;++br)
518 {
519 Update_PMat_At_Given_Edge(tree->a_edges[br],tree);
520 }
521
522 if(tree->n_root && tree->ignore_root == NO)
523 {
524 Update_PMat_At_Given_Edge(tree->n_root->b[1],tree);
525 Update_PMat_At_Given_Edge(tree->n_root->b[2],tree);
526 }
527 }
528 else//Update PMat for a specific edge
529 {
530 if(tree->use_eigen_lr == NO) Update_PMat_At_Given_Edge(b,tree);
531 }
532
533 if(!b)
534 {
535 if(tree->n_root != NULL)
536 {
537 if(tree->ignore_root == NO)
538 {
539 Post_Order_Lk(tree->n_root,tree->n_root->v[1],tree);
540 Post_Order_Lk(tree->n_root,tree->n_root->v[2],tree);
541
542 Update_Partial_Lk(tree,tree->n_root->b[1],tree->n_root);
543 Update_Partial_Lk(tree,tree->n_root->b[2],tree->n_root);
544
545 if(tree->both_sides == YES)
546 {
547 Pre_Order_Lk(tree->n_root,tree->n_root->v[2],tree);
548 Pre_Order_Lk(tree->n_root,tree->n_root->v[1],tree);
549 }
550 }
551 else
552 {
553 Post_Order_Lk(tree->e_root->rght,tree->e_root->left,tree);
554 Post_Order_Lk(tree->e_root->left,tree->e_root->rght,tree);
555
556 if(tree->both_sides == YES)
557 {
558 Pre_Order_Lk(tree->e_root->rght,tree->e_root->left,tree);
559 Pre_Order_Lk(tree->e_root->left,tree->e_root->rght,tree);
560 }
561 }
562 }
563 else
564 {
565 Post_Order_Lk(tree->a_nodes[tree->tip_root],tree->a_nodes[tree->tip_root]->v[0],tree);
566 if(tree->both_sides == YES)
567 Pre_Order_Lk(tree->a_nodes[tree->tip_root],tree->a_nodes[tree->tip_root]->v[0],tree);
568 }
569 }
570 }
571
572 if(!b)
573 {
574 if(tree->n_root)
575 {
576 if(tree->ignore_root == NO)
577 b = (tree->n_root->v[1]->tax == NO)?(tree->n_root->b[2]):(tree->n_root->b[1]);
578 else
579 b = tree->e_root;
580 }
581 else
582 b = tree->a_nodes[tree->tip_root]->b[0];
583 }
584
585 tree->c_lnL = .0;
586 tree->sum_min_sum_scale = .0;
587
588 #ifdef BEAGLE
589 calc_edgelks_beagle(b, tree);
590 #else
591
592
593 if(tree->update_eigen_lr == YES) Update_Eigen_Lr(b,tree);
594
595 if(tree->use_eigen_lr == YES)
596 {
597 for(catg=0;catg<ncatg;++catg)
598 {
599 len = MAX(0.0,b->l->v)*tree->mod->ras->gamma_rr->v[catg];
600 len *= tree->mod->br_len_mult->v;
601 if(tree->mixt_tree != NULL) len *= tree->mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number];
602 if(len < tree->mod->l_min) len = tree->mod->l_min;
603 else if(len > tree->mod->l_max) len = tree->mod->l_max;
604 for(state=0;state<ns;++state) expl[catg*ns+state] = (phydbl)pow(tree->mod->eigen->e_val[state],len);
605 }
606 }
607
608 p_lk_left = b->p_lk_left;
609 p_lk_rght = b->rght->tax ? b->p_lk_tip_r : b->p_lk_rght;
610
611 CALL++;
612
613
614 for(site=0;site<npatterns;++site)
615 {
616 ambiguity_check = -1;
617 state = -1;
618 tree->curr_site = site;
619
620 if((b->rght->tax) && (tree->mod->s_opt->greedy == NO))
621 {
622 ambiguity_check = b->rght->c_seq->is_ambigu[tree->curr_site];
623 if(ambiguity_check == NO)
624 {
625 state = b->rght->c_seq->d_state[tree->curr_site];
626 }
627 }
628
629 if(tree->mod->use_m4mod) ambiguity_check = YES;
630
631 if(tree->use_eigen_lr == YES)
632 {
633 if(tree->data->wght[site] > SMALL) Lk_Core_Eigen_Lr(expl,dot_prod,b,tree);
634 dot_prod += nsncatg;
635 }
636 else
637 {
638 if(tree->data->wght[site] > SMALL) Lk_Core(state,ambiguity_check,p_lk_left,p_lk_rght,b->Pij_rr,b->tPij_rr,b,tree);
639
640 if(b->rght->tax == YES)
641 {
642 p_lk_left += nsncatg;
643 p_lk_rght += ns;
644 }
645 else
646 {
647 p_lk_left += nsncatg;
648 p_lk_rght += nsncatg;
649 }
650 }
651 }
652 #endif
653
654 return tree->c_lnL;
655 }
656
657 //////////////////////////////////////////////////////////////
658 //////////////////////////////////////////////////////////////
659 // First derivative of the log-likelihood with respect
660 // to the length of edge b
dLk(phydbl * l,t_edge * b,t_tree * tree)661 phydbl dLk(phydbl *l, t_edge *b, t_tree *tree)
662 {
663 unsigned int catg,state,site;
664 phydbl len,rr;
665 phydbl lk,dlk,dlnlk,lnlk;
666 phydbl ev,expevlen;
667
668 const unsigned int ns = tree->mod->ns;
669 const unsigned int ncatg = tree->mod->ras->n_catg;
670 const unsigned int npattern = tree->n_pattern;
671
672 phydbl *dot_prod = tree->dot_prod;
673 phydbl *expl = tree->expl;
674
675 tree->numerical_warning = NO;
676
677 assert(isnan(*l) == FALSE);
678
679 if(*l < tree->mod->l_min) *l = tree->mod->l_min;
680 else if(*l > tree->mod->l_max) *l = tree->mod->l_max;
681
682 assert(b != NULL);
683
684 if(tree->is_mixt_tree == YES)
685 {
686 #ifdef BEAGLE
687 Warn_And_Exit(TODO_BEAGLE);
688 #endif
689 return MIXT_dLk(l,b,tree);
690 }
691
692 if(tree->update_eigen_lr == YES) Update_Eigen_Lr(b,tree);
693
694 for(catg=0;catg<ncatg;catg++)
695 {
696 rr = tree->mod->ras->gamma_rr->v[catg];
697 rr *= tree->mod->br_len_mult->v;
698
699 len = (*l) * rr;
700
701 if(isinf(len) || isnan(len))
702 {
703 PhyML_Fprintf(stderr,"\n. len=%f rr=%f l=%f",len,rr,*l);
704 assert(FALSE);
705 }
706
707 if(len < tree->mod->l_min) len = tree->mod->l_min;
708 else if(len > tree->mod->l_max) len = tree->mod->l_max;
709 // value of rr should be corrected too if any of these two conditions
710 // is true. Leads to numerical precision issues though...
711
712
713 for(state=0;state<ns;++state)
714 {
715 ev = tree->mod->eigen->e_val[state];
716 expevlen = exp(ev*len);
717
718 expl[catg*2*ns + 2*state] = expevlen;
719 expl[catg*2*ns + 2*state + 1] = expevlen*ev*rr;
720 }
721 }
722
723 dlnlk = 0.0;
724 lnlk = 0.0;
725
726 for(site=0;site<npattern;++site)
727 {
728 if(tree->data->wght[site] > SMALL)
729 {
730 tree->curr_site = site;
731
732 Lk_dLk_Core_Eigen_Lr(expl,dot_prod+site*ns*ncatg,b,&lk,&dlk,tree);
733
734 assert(lk > .0);
735
736 dlk /= lk;
737 dlnlk += tree->data->wght[site] * dlk;
738 lnlk += tree->data->wght[site] * (log(lk) - (phydbl)LOG2 * tree->fact_sum_scale[site]);
739 }
740 }
741
742 tree->c_dlnL = dlnlk;
743 tree->c_lnL = lnlk;
744
745 return tree->c_lnL;
746 }
747
748 //////////////////////////////////////////////////////////////
749 //////////////////////////////////////////////////////////////
750
751 /* Core of the likelihood calculation. Assume that the partial likelihoods on both
752 sides of t_edge *b are up-to-date. Calculate the log-likelihood at one site.
753 Note: this function can be used to evaluate first or second derivative of the
754 likelihood function with respect to the length of b, at a given site. Hence, be
755 careful with the meaning of 'site_lk'. If 'derivative=TRUE', then site_lk is
756 either the first or second derivative of the likelihood at that site, given the
757 length of edge 'b'.
758 */
759
Lk_Core(int state,int ambiguity_check,phydbl * p_lk_left,phydbl * p_lk_rght,phydbl * Pij_rr,phydbl * tPij_rr,t_edge * b,t_tree * tree)760 phydbl Lk_Core(int state, int ambiguity_check,
761 phydbl *p_lk_left, phydbl *p_lk_rght,
762 phydbl *Pij_rr,
763 phydbl *tPij_rr,
764 t_edge *b,
765 t_tree *tree)
766 {
767 phydbl site_lk,res,*pi,*site_lk_cat,log_site_lk;
768 unsigned int catg;
769
770 const unsigned int ns = tree->mod->ns;
771 const unsigned int ncatg = tree->mod->ras->n_catg;
772 const unsigned int site = tree->curr_site;
773 const unsigned nsns = ns*ns;
774
775
776 assert(tree->data->wght[site] > SMALL);
777
778 pi = tree->mod->e_frq->pi->v;
779
780 if(tree->mod->s_opt->skip_tree_traversal == NO)
781 {
782 for(catg=0;catg<ncatg;++catg)
783 {
784 if(ns == 4 || ns == 20)
785 {
786 #if (defined(__AVX__) || defined(__AVX2__))
787 tree->site_lk_cat[catg] = AVX_Lk_Core_One_Class_No_Eigen_Lr(p_lk_left,p_lk_rght,Pij_rr,tPij_rr,pi,ns,ambiguity_check,state);
788 #elif (defined(__SSE3__))
789 tree->site_lk_cat[catg] = SSE_Lk_Core_One_Class_No_Eigen_Lr(p_lk_left,p_lk_rght,Pij_rr,tPij_rr,pi,ns,ambiguity_check,state);
790 #else
791 tree->site_lk_cat[catg] = Lk_Core_One_Class_No_Eigen_Lr(p_lk_left,p_lk_rght,Pij_rr,pi,ns,ambiguity_check,state);
792 #endif
793 }
794 else
795 {
796 tree->site_lk_cat[catg] = Lk_Core_One_Class_No_Eigen_Lr(p_lk_left,p_lk_rght,Pij_rr,pi,ns, ambiguity_check, state);
797 }
798
799 Pij_rr += nsns;
800 tPij_rr += nsns;
801 if(b->left->tax == NO) p_lk_left += ns;
802 if(b->rght->tax == NO) p_lk_rght += ns;
803 }
804
805 Pull_Scaling_Factors(site,b,tree);
806
807 }
808
809 site_lk = .0;
810 site_lk_cat = tree->unscaled_site_lk_cat + site*ncatg;
811 for(catg=0;catg<ncatg;++catg) site_lk += site_lk_cat[catg] * tree->mod->ras->gamma_r_proba->v[catg];
812
813 if(tree->mod->ras->invar == YES)
814 {
815 int num_prec_issue = NO;
816 phydbl inv_site_lk = Invariant_Lk(tree->fact_sum_scale[site],site,&num_prec_issue,tree);
817
818 switch(num_prec_issue)
819 {
820 case YES :
821 {
822 assert(isinf(inv_site_lk));
823 tree->fact_sum_scale[site] = 0;
824 inv_site_lk = Invariant_Lk(0,site,&num_prec_issue,tree);
825 site_lk = inv_site_lk * tree->mod->ras->pinvar->v;
826 break;
827 }
828 case NO :
829 {
830 site_lk = site_lk * (1. - tree->mod->ras->pinvar->v) + inv_site_lk * tree->mod->ras->pinvar->v;
831 break;
832 }
833 }
834
835 }
836
837 if(tree->apply_lk_scaling == YES) res = site_lk / pow(2,tree->fact_sum_scale[site]);
838 else res = site_lk;
839
840 if(site_lk < SMALL)
841 {
842 site_lk = SMALL;
843 tree->numerical_warning = YES;
844 }
845
846
847 log_site_lk = log(site_lk) - (phydbl)LOG2 * tree->fact_sum_scale[site]; // log_site_lk = log(site_lk_scaled / 2^(left_subtree+right_subtree))
848 tree->c_lnL_sorted[site] = log_site_lk;
849 tree->c_lnL += tree->data->wght[site] * log_site_lk;
850 tree->cur_site_lk[site] = exp(log_site_lk); // note to self : add opt out option to avoid calculating this if not necessary
851
852 /* printf("\n. clnL: %f",tree->c_lnL); */
853 return res;
854 }
855
856 //////////////////////////////////////////////////////////////
857 //////////////////////////////////////////////////////////////
858
Lk_Core_Eigen_Lr(phydbl * expl,phydbl * dot_prod,t_edge * b,t_tree * tree)859 phydbl Lk_Core_Eigen_Lr(phydbl *expl, phydbl *dot_prod, t_edge *b, t_tree *tree)
860 {
861 phydbl site_lk,res,*site_lk_cat,log_site_lk;
862 unsigned int catg;
863 int num_prec_issue;
864
865 const unsigned int ns = tree->mod->ns;
866 const unsigned int ncatg = tree->mod->ras->n_catg;
867 const unsigned int site = tree->curr_site;
868
869 assert(tree->data->wght[site] > SMALL);
870
871 if(tree->mod->s_opt->skip_tree_traversal == NO)
872 {
873 for(catg=0;catg<ncatg;++catg)
874 {
875 if(tree->mod->io->datatype == NT || tree->mod->io->datatype == AA)
876 {
877 #if (defined(__AVX__) || defined(__AVX2__))
878 tree->site_lk_cat[catg] = AVX_Lk_Core_One_Class_Eigen_Lr(dot_prod,expl ? expl : NULL,ns);
879 #elif (defined(__SSE3__))
880 tree->site_lk_cat[catg] = SSE_Lk_Core_One_Class_Eigen_Lr(dot_prod,expl ? expl : NULL,ns);
881 #else
882 tree->site_lk_cat[catg] = Lk_Core_One_Class_Eigen_Lr(dot_prod,expl ? expl : NULL,ns);
883 #endif
884 }
885 else
886 {
887 tree->site_lk_cat[catg] = Lk_Core_One_Class_Eigen_Lr(dot_prod,expl ? expl : NULL,ns);
888 }
889
890 dot_prod += ns;
891 if(expl) expl += ns;
892 }
893
894 Pull_Scaling_Factors(site,b,tree);
895
896 }
897
898
899 site_lk_cat = tree->unscaled_site_lk_cat + site*ncatg;
900 site_lk = .0;
901 for(catg=0;catg<ncatg;++catg) site_lk += site_lk_cat[catg] * tree->mod->ras->gamma_r_proba->v[catg];
902
903 if(tree->mod->ras->invar == YES)
904 {
905 num_prec_issue = NO;
906 phydbl inv_site_lk = Invariant_Lk(tree->fact_sum_scale[site],site,&num_prec_issue,tree);
907
908 switch(num_prec_issue)
909 {
910 case YES :
911 {
912 assert(isinf(inv_site_lk));
913 tree->fact_sum_scale[site] = 0;
914 inv_site_lk = Invariant_Lk(0,site,&num_prec_issue,tree);
915 site_lk = inv_site_lk * tree->mod->ras->pinvar->v;
916 break;
917 }
918 case NO :
919 {
920 site_lk = site_lk * (1. - tree->mod->ras->pinvar->v) + inv_site_lk * tree->mod->ras->pinvar->v;
921 break;
922 }
923 }
924 }
925
926 if(site_lk < SMALL)
927 {
928 site_lk = SMALL;
929 tree->numerical_warning = YES;
930 }
931
932 // likelihood (or 1st, 2nd derivative) not rescaled here. Valid only if all partial likelihoods
933 // were scaled using the same factor, i.e., when scaling_method == SCALE_FAST. In this case, the
934 // scaling factors will cancel out in dlk/lk and d2lk/lk
935 res = site_lk;
936
937 log_site_lk = log(site_lk) - (phydbl)LOG2 * tree->fact_sum_scale[site]; // log_site_lk = log(site_lk_scaled / 2^(left_subtree+right_subtree))
938 tree->c_lnL_sorted[site] = log_site_lk;
939 tree->c_lnL += tree->data->wght[site] * log_site_lk;
940 tree->cur_site_lk[site] = exp(log_site_lk); // note to self : add opt out option to avoid calculating this if not necessary
941
942 return res;
943 }
944
945 //////////////////////////////////////////////////////////////
946 //////////////////////////////////////////////////////////////
947 // Compute likelihood and first derivative of likelihood with respect to the length of edge b *unscaled*
Lk_dLk_Core_Eigen_Lr(phydbl * expl,phydbl * dot_prod,t_edge * b,phydbl * lk,phydbl * dlk,t_tree * tree)948 void Lk_dLk_Core_Eigen_Lr(phydbl *expl, phydbl *dot_prod, t_edge *b, phydbl *lk, phydbl *dlk, t_tree *tree)
949 {
950 phydbl core_lk,core_dlk;
951 unsigned int catg;
952 int num_prec_issue;
953
954 const unsigned int ns = tree->mod->ns;
955 const unsigned int ncatg = tree->mod->ras->n_catg;
956 const unsigned int site = tree->curr_site;
957
958 *lk = *dlk = 0.0;
959
960 assert(tree->data->wght[site] > SMALL);
961
962 if(tree->mod->s_opt->skip_tree_traversal == NO)
963 {
964 for(catg=0;catg<ncatg;++catg)
965 {
966 if(tree->mod->io->datatype == NT || tree->mod->io->datatype == AA)
967 {
968 #if (defined(__AVX__) || defined(__AVX2__))
969 AVX_Lk_dLk_Core_One_Class_Eigen_Lr(dot_prod,
970 expl ? expl : NULL,
971 ns,&core_lk,&core_dlk);
972 #elif (defined(__SSE3__))
973 SSE_Lk_dLk_Core_One_Class_Eigen_Lr(dot_prod,
974 expl ? expl : NULL,
975 ns,&core_lk,&core_dlk);
976 #else
977 Lk_dLk_Core_One_Class_Eigen_Lr(dot_prod,
978 expl ? expl : NULL,
979 ns,&core_lk,&core_dlk);
980 #endif
981 }
982 else
983 {
984 Lk_dLk_Core_One_Class_Eigen_Lr(dot_prod,
985 expl ? expl : NULL,
986 ns,&core_lk,&core_dlk);
987 }
988
989 *lk += core_lk * tree->mod->ras->gamma_r_proba->v[catg];
990 *dlk += core_dlk * tree->mod->ras->gamma_r_proba->v[catg];
991
992 dot_prod += ns;
993 if(expl) expl += 2*ns;
994 }
995 Pull_Scaling_Factors(site,b,tree);
996 }
997
998 if(tree->mod->ras->invar == YES)
999 {
1000 num_prec_issue = NO;
1001 phydbl inv_site_lk = Invariant_Lk(tree->fact_sum_scale[site],site,&num_prec_issue,tree);
1002
1003 switch(num_prec_issue)
1004 {
1005 case YES :
1006 {
1007 *lk = inv_site_lk * tree->mod->ras->pinvar->v;
1008 *dlk = 0.0;
1009 break;
1010 }
1011 case NO :
1012 {
1013 *lk = *lk * (1. - tree->mod->ras->pinvar->v) + inv_site_lk * tree->mod->ras->pinvar->v;
1014 *dlk = *dlk * (1. - tree->mod->ras->pinvar->v);
1015 break;
1016 }
1017 }
1018 }
1019
1020 if(*lk < SMALL)
1021 {
1022 *lk = SMALL;
1023 tree->numerical_warning = YES;
1024 }
1025 }
1026
1027
1028 //////////////////////////////////////////////////////////////
1029 //////////////////////////////////////////////////////////////
1030
Update_Eigen_Lr(t_edge * b,t_tree * tree)1031 void Update_Eigen_Lr(t_edge *b, t_tree *tree)
1032 {
1033 unsigned int site,catg,i,j;
1034 phydbl *dot_prod,*r_e_vect,*l_e_vect,*p_lk_left,*p_lk_rght,*pi;
1035 phydbl left,rght;
1036
1037 const unsigned int npattern = tree->n_pattern;
1038 const unsigned int ns = tree->mod->ns;
1039 const unsigned int ncatg = tree->mod->ras->n_catg;
1040 const unsigned int nsncatg = ns*ncatg;
1041
1042
1043 if(tree->is_mixt_tree == YES)
1044 {
1045 MIXT_Update_Eigen_Lr(b,tree);
1046 return;
1047 }
1048
1049 if(tree->mod->ns == 4 || tree->mod->ns == 20)
1050 {
1051 #if (defined(__AVX__) || defined(__AVX2__))
1052 AVX_Update_Eigen_Lr(b,tree);
1053 return;
1054 #elif (defined(__SSE3__))
1055 SSE_Update_Eigen_Lr(b,tree);
1056 return;
1057 #endif
1058 }
1059
1060 assert(tree->update_eigen_lr == YES);
1061
1062 dot_prod = tree->dot_prod;
1063 r_e_vect = tree->mod->eigen->r_e_vect;
1064 l_e_vect = tree->mod->eigen->l_e_vect;
1065 pi = tree->mod->e_frq->pi->v;
1066
1067 if(b->left->tax == YES) p_lk_left = b->p_lk_tip_l;
1068 else p_lk_left = b->p_lk_left;
1069
1070 if(b->rght->tax == YES) p_lk_rght = b->p_lk_tip_r;
1071 else p_lk_rght = b->p_lk_rght;
1072
1073 for(site=0;site<npattern;++site)
1074 {
1075 if(tree->data->wght[site] > SMALL)
1076 {
1077 for(catg=0;catg<ncatg;++catg)
1078 {
1079 for(i=0;i<ns;++i)
1080 {
1081 left = rght = 0.0;
1082 for(j=0;j<ns;++j)
1083 {
1084 left += r_e_vect[j*ns + i] * p_lk_left[j] * pi[j];
1085 rght += l_e_vect[i*ns + j] * p_lk_rght[j];
1086 }
1087 dot_prod[i] = left*rght;
1088 }
1089 dot_prod += ns;
1090 if(b->left->tax == NO) p_lk_left += ns;
1091 if(b->rght->tax == NO) p_lk_rght += ns;
1092 }
1093 if(b->left->tax == YES) p_lk_left += ns;
1094 if(b->rght->tax == YES) p_lk_rght += ns;
1095 }
1096 else
1097 {
1098 if(b->left->tax == YES) p_lk_left += ns;
1099 else p_lk_left += nsncatg;
1100
1101 if(b->rght->tax == YES) p_lk_rght += ns;
1102 else p_lk_rght += nsncatg;
1103
1104 dot_prod += nsncatg;
1105 }
1106 }
1107 }
1108
1109 //////////////////////////////////////////////////////////////
1110 //////////////////////////////////////////////////////////////
1111
Rate_Correction(int exponent,phydbl * site_lk_cat)1112 void Rate_Correction(int exponent, phydbl *site_lk_cat)
1113 {
1114 int piecewise_exponent;
1115 phydbl multiplier,dum;
1116 unsigned long long int one = 1;
1117
1118 dum = *site_lk_cat;
1119 if(exponent >= 0)
1120 {
1121 /* Multiply by 2^exponent */
1122 do
1123 {
1124 piecewise_exponent = MIN(exponent,63);
1125 multiplier = (phydbl)(one << piecewise_exponent);
1126 dum = dum * multiplier;
1127 exponent = exponent - piecewise_exponent;
1128 }
1129 while(exponent != 0);
1130 }
1131 else
1132 {
1133 /* Divide by 2^exponent */
1134 do
1135 {
1136 piecewise_exponent = MAX(exponent,-63);
1137 multiplier = 1. / (phydbl)(one << -piecewise_exponent);
1138 dum = dum * multiplier;
1139 exponent = exponent - piecewise_exponent;
1140 }
1141 while(exponent != 0);
1142 }
1143
1144 *site_lk_cat = dum;
1145 }
1146
1147 //////////////////////////////////////////////////////////////
1148 //////////////////////////////////////////////////////////////
1149
Lk_Core_One_Class_Eigen_Lr(phydbl * dot_prod,phydbl * expl,int ns)1150 phydbl Lk_Core_One_Class_Eigen_Lr(phydbl *dot_prod, phydbl *expl, int ns)
1151 {
1152 unsigned int l;
1153 phydbl lk = 0.0;
1154 if(expl != NULL) for(l=0;l<ns;++l) lk += dot_prod[l] * expl[l];
1155 else for(l=0;l<ns;++l) lk += dot_prod[l];
1156 return lk;
1157 }
1158
1159
1160 //////////////////////////////////////////////////////////////
1161 //////////////////////////////////////////////////////////////
1162
Lk_dLk_Core_One_Class_Eigen_Lr(phydbl * dot_prod,phydbl * expl,unsigned int ns,phydbl * lk,phydbl * dlk)1163 void Lk_dLk_Core_One_Class_Eigen_Lr(phydbl *dot_prod, phydbl *expl, unsigned int ns, phydbl *lk, phydbl *dlk)
1164 {
1165 unsigned int i;
1166
1167 *lk = *dlk = 0.0;
1168 for(i=0;i<ns;++i)
1169 {
1170 *lk += dot_prod[i] * expl[2*i];
1171 *dlk += dot_prod[i] * expl[2*i+1];
1172 }
1173 }
1174
1175 //////////////////////////////////////////////////////////////
1176 //////////////////////////////////////////////////////////////
1177
Lk_Core_One_Class_No_Eigen_Lr(phydbl * p_lk_left,phydbl * p_lk_rght,phydbl * Pij,phydbl * pi,int ns,int ambiguity_check,int state)1178 phydbl Lk_Core_One_Class_No_Eigen_Lr(phydbl *p_lk_left, phydbl *p_lk_rght, phydbl *Pij,phydbl *pi, int ns, int ambiguity_check, int state)
1179 {
1180 unsigned int l,k;
1181 phydbl lk = 0.0;
1182 phydbl sum;
1183
1184
1185 if(ambiguity_check == NO)/* tip case */
1186 {
1187 sum = .0;
1188 Pij += state*ns;
1189 for(l=0;l<ns;++l) sum += Pij[l] * p_lk_left[l];
1190 lk += sum * pi[state];
1191 }
1192 else /* If the character observed at the tip is ambiguous: ns x ns terms to consider */
1193 {
1194 for(k=0;k<ns;++k)
1195 {
1196 if(p_lk_rght[k] > .0) /* Only bother ascending into the subtrees if the likelihood of state k, at site "site*dim2" is > 0 */
1197 {
1198 sum = .0;
1199 for(l=0;l<ns;l++)
1200 {
1201 sum += Pij[l] * p_lk_left[l];
1202 }
1203 lk += sum * pi[k] * p_lk_rght[k];
1204 }
1205 Pij += ns;
1206 }
1207 }
1208
1209 return lk;
1210 }
1211
1212 /* #endif */
1213
1214 //////////////////////////////////////////////////////////////
1215 //////////////////////////////////////////////////////////////
1216
1217 // Returns the scaled likelihood for invariable sites
Invariant_Lk(int fact_sum_scale,int site,int * num_prec_issue,t_tree * tree)1218 phydbl Invariant_Lk(int fact_sum_scale, int site, int *num_prec_issue, t_tree *tree)
1219 {
1220 int exponent,piecewise_exponent;
1221 phydbl multiplier;
1222 phydbl inv_site_lk = 0.;
1223
1224 (*num_prec_issue) = NO;
1225
1226 /* The substitution model does include invariable sites */
1227 if(tree->mod->ras->invar == YES)
1228 {
1229 /* The site is invariant */
1230 if(tree->data->invar[site] > -0.5)
1231 {
1232 inv_site_lk = tree->mod->e_frq->pi->v[tree->data->invar[site]];
1233
1234 /* printf("\n. inv_site_lk = %f [%c] [%d] invar: %d",inv_site_lk,tree->data->c_seq[0]->state[site],tree->data->invar[site],tree->data->invar[site]); */
1235
1236 if(tree->apply_lk_scaling == YES)
1237 {
1238 exponent = fact_sum_scale;
1239 do
1240 {
1241 piecewise_exponent = MIN(exponent,63);
1242 multiplier = (phydbl)((unsigned long long)(1) << piecewise_exponent);
1243 inv_site_lk *= multiplier;
1244 exponent -= piecewise_exponent;
1245 }
1246 while(exponent != 0);
1247 }
1248
1249 /* Update the value of site_lk */
1250 if(isinf(inv_site_lk)) // P(D|r=0) >> P(D|r>0) => assume P(D) = P(D|r=0)P(r=0)
1251 {
1252 int i;
1253 PhyML_Fprintf(stderr,"\n. fact_sum_scale: %d",fact_sum_scale);
1254 PhyML_Fprintf(stderr,"\n. pi: %f",tree->mod->e_frq->pi->v[tree->data->invar[site]]);
1255 for(i=0;i<tree->mod->ns;i++) PhyML_Fprintf(stderr,"\n. pi %d: %f",i,tree->mod->e_frq->pi->v[i]);
1256 PhyML_Fprintf(stderr,"\n. Numerical precision issue alert.");
1257 PhyML_Fprintf(stderr,"\n. File %s at line %d (function '%s')\n",__FILE__,__LINE__,__FUNCTION__);
1258 (*num_prec_issue) = YES;
1259 }
1260 }
1261 }
1262
1263 return inv_site_lk;
1264
1265 }
1266
1267 //////////////////////////////////////////////////////////////
1268 //////////////////////////////////////////////////////////////
1269
1270 /* Update partial likelihood on edge b on the side of b where
1271 node d lies.
1272 */
1273
Update_Partial_Lk(t_tree * tree,t_edge * b,t_node * d)1274 void Update_Partial_Lk(t_tree *tree, t_edge *b, t_node *d)
1275 {
1276 /* if(tree->eval_alnL == NO) return; */
1277 if(b->left == d && b->update_partial_lk_left == NO) return;
1278 if(b->rght == d && b->update_partial_lk_rght == NO) return;
1279
1280 if(tree->is_mixt_tree)
1281 {
1282 MIXT_Update_Partial_Lk(tree,b,d);
1283 return;
1284 }
1285
1286 if((tree->io->do_alias_subpatt == YES) &&
1287 (tree->update_alias_subpatt == YES))
1288 Alias_One_Subpatt((d==b->left)?(b->rght):(b->left),d,tree);
1289
1290 if(d->tax) return;
1291
1292
1293 #ifdef BEAGLE
1294 update_beagle_partials(tree, b, d);
1295 #else
1296 if(tree->mod->use_m4mod == NO)
1297 {
1298 if(tree->mod->ns == 4 || tree->mod->ns == 20)
1299 {
1300 #if (defined(__AVX__) || defined(__AVX2__))
1301 AVX_Update_Partial_Lk(tree,b,d);
1302 #elif (defined(__SSE3__))
1303 SSE_Update_Partial_Lk(tree,b,d);
1304 #else
1305 Default_Update_Partial_Lk(tree,b,d);
1306 #endif
1307 }
1308 else
1309 {
1310 Update_Partial_Lk_Generic(tree,b,d);
1311 }
1312 }
1313 else
1314 {
1315 Update_Partial_Lk_Generic(tree,b,d);
1316 }
1317 #endif
1318
1319
1320 // Print_Edge_Likelihoods(tree, b, false);
1321 }
1322
1323 //////////////////////////////////////////////////////////////
1324 //////////////////////////////////////////////////////////////
1325
1326 #ifndef BEAGLE
1327
Update_Partial_Lk_Generic(t_tree * tree,t_edge * b,t_node * d)1328 void Update_Partial_Lk_Generic(t_tree *tree, t_edge *b, t_node *d)
1329 {
1330 /*
1331 |
1332 |<- b
1333 |
1334 d
1335 / \
1336 / \
1337 / \
1338 n_v1 n_v2
1339 */
1340 t_node *n_v1, *n_v2;
1341 phydbl p1_lk1,p2_lk2;
1342 phydbl *p_lk,*p_lk_v1,*p_lk_v2;
1343 phydbl *Pij1,*Pij2;
1344 phydbl *tPij1,*tPij2;
1345 int *sum_scale, *sum_scale_v1, *sum_scale_v2;
1346 int sum_scale_v1_val, sum_scale_v2_val;
1347 int i,j;
1348 int catg,site;
1349 int n_patterns;
1350 short int ambiguity_check_v1,ambiguity_check_v2;
1351 int state_v1,state_v2;
1352 phydbl smallest_p_lk,largest_p_lk;
1353 int *p_lk_loc;
1354
1355 unsigned const int ncatg = tree->mod->ras->n_catg;
1356 unsigned const int ns = tree->mod->ns;
1357 unsigned const int ncatgns = ncatg * ns;
1358 unsigned const int nsns = ns * ns;
1359
1360
1361
1362 if(tree->n_root && tree->ignore_root == YES &&
1363 (d == tree->n_root->v[1] || d == tree->n_root->v[2]) &&
1364 (b == tree->n_root->b[1] || b == tree->n_root->b[2]))
1365 {
1366 assert(FALSE);
1367 }
1368
1369
1370 state_v1 = state_v2 = -1;
1371 ambiguity_check_v1 = ambiguity_check_v2 = NO;
1372 sum_scale_v1_val = sum_scale_v2_val = 0;
1373 p1_lk1 = p2_lk2 = .0;
1374
1375 if(d->tax)
1376 {
1377 PhyML_Fprintf(stderr,"\n. t_node %d is a leaf...",d->num);
1378 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s').\n",__FILE__,__LINE__,__FUNCTION__);
1379 Exit("\n");
1380 }
1381
1382 n_patterns = tree->n_pattern;
1383
1384 n_v1 = n_v2 = NULL;
1385 p_lk = p_lk_v1 = p_lk_v2 = NULL;
1386 Pij1 = Pij2 = NULL;
1387 tPij1 = tPij2 = NULL;
1388 sum_scale_v1 = sum_scale_v2 = NULL;
1389 p_lk_loc = NULL;
1390 smallest_p_lk = BIG;
1391
1392 Set_All_Partial_Lk(&n_v1,&n_v2,
1393 &p_lk,&sum_scale,&p_lk_loc,
1394 &Pij1,&tPij1,&p_lk_v1,&sum_scale_v1,
1395 &Pij2,&tPij2,&p_lk_v2,&sum_scale_v2,
1396 d,b,tree);
1397
1398 /* For every site in the alignment */
1399 for(site=0;site<n_patterns;site++)
1400 {
1401 if(tree->data->wght[site] > SMALL)
1402 {
1403 state_v1 = state_v2 = -1;
1404 ambiguity_check_v1 = ambiguity_check_v2 = NO;
1405
1406 if(tree->mod->s_opt->greedy == NO)
1407 {
1408 /* n_v1 and n_v2 are tip nodes */
1409 if(n_v1 && n_v1->tax)
1410 {
1411 /* Is the state at this tip ambiguous? */
1412 ambiguity_check_v1 = n_v1->c_seq->is_ambigu[site];
1413 /* if(ambiguity_check_v1 == NO) state_v1 = Get_State_From_Partial_Pars(n_v1->b[0]->p_lk_tip_r,site*ns,tree); */
1414 if(ambiguity_check_v1 == NO) state_v1 = n_v1->c_seq->d_state[site];
1415 }
1416
1417 if(n_v2 && n_v2->tax)
1418 {
1419 /* Is the state at this tip ambiguous? */
1420 ambiguity_check_v2 = n_v2->c_seq->is_ambigu[site];
1421 /* ambiguity_check_v2 = tree->data->c_seq[n_v2->num]->is_ambigu[site]; */
1422 /* if(ambiguity_check_v2 == NO) state_v2 = Get_State_From_Partial_Pars(n_v2->b[0]->p_lk_tip_r,site*ns,tree); */
1423 if(ambiguity_check_v2 == NO) state_v2 = n_v2->c_seq->d_state[site];
1424 }
1425 }
1426
1427 if(tree->mod->use_m4mod)
1428 {
1429 ambiguity_check_v1 = YES;
1430 ambiguity_check_v2 = YES;
1431 }
1432
1433 /* For all the rate classes */
1434 for(catg=0;catg<ncatg;catg++)
1435 {
1436 if(tree->mod->ras->skip_rate_cat[catg] == YES) continue;
1437
1438 smallest_p_lk = BIG;
1439
1440 /* For all the states at node d */
1441 for(i=0;i<tree->mod->ns;i++)
1442 {
1443 p1_lk1 = .0;
1444
1445 if(n_v1)
1446 {
1447 /* n_v1 is a tip */
1448 if((n_v1->tax) && (!tree->mod->s_opt->greedy))
1449 {
1450 if(ambiguity_check_v1 == NO)
1451 {
1452 /* For the (non-ambiguous) state at node n_v1 */
1453 p1_lk1 = Pij1[catg*nsns+i*ns+state_v1];
1454 }
1455 else
1456 {
1457 /* For all the states at node n_v1 */
1458 for(j=0;j<tree->mod->ns;j++)
1459 {
1460 p1_lk1 += Pij1[catg*nsns+i*ns+j] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*ns+j];
1461 }
1462 }
1463 }
1464 /* n_v1 is an internal node */
1465 else
1466 {
1467 /* For the states at node n_v1 */
1468 for(j=0;j<tree->mod->ns;j++)
1469 {
1470 p1_lk1 += Pij1[catg*nsns+i*ns+j] * p_lk_v1[site*ncatgns+catg*ns+j];
1471 }
1472 }
1473 }
1474 else
1475 {
1476 p1_lk1 = 1.0;
1477 }
1478
1479 p2_lk2 = .0;
1480
1481 /* We do exactly the same as for node n_v1 but for node n_v2 this time.*/
1482 if(n_v2)
1483 {
1484 /* n_v2 is a tip */
1485 if((n_v2->tax) && (!tree->mod->s_opt->greedy))
1486 {
1487 if(ambiguity_check_v2 == NO)
1488 {
1489 /* For the (non-ambiguous) state at node n_v2 */
1490 p2_lk2 = Pij2[catg*nsns+i*ns+state_v2];
1491 }
1492 else
1493 {
1494 /* For all the states at node n_v2 */
1495 for(j=0;j<tree->mod->ns;j++)
1496 {
1497 p2_lk2 += Pij2[catg*nsns+i*ns+j] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*ns+j];
1498 }
1499 }
1500 }
1501 /* n_v2 is an internal node */
1502 else
1503 {
1504 /* For all the states at node n_v2 */
1505 for(j=0;j<tree->mod->ns;j++)
1506 {
1507 p2_lk2 += Pij2[catg*nsns+i*ns+j] * p_lk_v2[site*ncatgns+catg*ns+j];
1508 }
1509 }
1510 }
1511 else
1512 {
1513 p2_lk2 = 1.0;
1514 }
1515
1516 p_lk[site*ncatgns+catg*ns+i] = p1_lk1 * p2_lk2;
1517
1518 /* if(site == 0) PhyML_Printf("\n+ site: %d %G",site,p_lk[site*ncatgns+catg*ns+i]); */
1519
1520 }
1521
1522 if(tree->scaling_method == SCALE_RATE_SPECIFIC)
1523 {
1524 smallest_p_lk = BIG;
1525 for(i=0;i<ns;++i)
1526 if(p_lk[site*ncatgns+catg*ns+i] < smallest_p_lk)
1527 smallest_p_lk = p_lk[site*ncatgns+catg*ns+i];
1528
1529 /* Current scaling values at that site */
1530 sum_scale_v1_val = (sum_scale_v1)?(sum_scale_v1[site*ncatg+catg]):(0);
1531 sum_scale_v2_val = (sum_scale_v2)?(sum_scale_v2[site*ncatg+catg]):(0);
1532
1533 sum_scale[site*ncatg+catg] = sum_scale_v1_val + sum_scale_v2_val;
1534
1535 /* Scaling. We have p_lk_lim_inf = 2^-500. Consider for instance that
1536 smallest_p_lk = 2^-600, then curr_scaler_pow will be equal to 100, and
1537 each element in the partial likelihood vector will be multiplied by
1538 2^100. */
1539 if(smallest_p_lk < (phydbl)P_LK_LIM_INF &&
1540 tree->mod->augmented == NO &&
1541 tree->apply_lk_scaling == YES &&
1542 (n_v1->tax == NO || n_v2->tax == NO))
1543
1544 {
1545 int curr_scaler_pow;
1546 curr_scaler_pow = (int)(-500.*LOG2-log(smallest_p_lk))/LOG2;
1547 sum_scale[site*ncatg+catg] += curr_scaler_pow;
1548 for(i=0;i<ns;++i) Rate_Correction(curr_scaler_pow, p_lk + site*nsns + catg*ns + i);
1549
1550 }
1551 }
1552 }
1553
1554 if(tree->scaling_method == SCALE_FAST)
1555 {
1556 sum_scale_v1_val = (sum_scale_v1)?(sum_scale_v1[site]):(0);
1557 sum_scale_v2_val = (sum_scale_v2)?(sum_scale_v2[site]):(0);
1558 sum_scale[site] = sum_scale_v1_val + sum_scale_v2_val;
1559
1560 sum_scale[site*ncatg+catg] = sum_scale_v1_val + sum_scale_v2_val;
1561
1562 assert(sum_scale[site] < 1024);
1563
1564 largest_p_lk = -BIG;
1565 for(i=0;i<ns*ncatg;++i)
1566 if(p_lk[site*ncatgns+i] > largest_p_lk)
1567 largest_p_lk = p_lk[site*ncatgns+i] ;
1568
1569 if(largest_p_lk < INV_TWO_TO_THE_LARGE &&
1570 tree->mod->augmented == NO &&
1571 tree->apply_lk_scaling == YES)
1572 {
1573 for(i=0;i<ns*ncatg;++i) p_lk[site*ncatgns + i] *= TWO_TO_THE_LARGE;
1574 sum_scale[site] += LARGE;
1575 }
1576 }
1577 }
1578 else
1579 {
1580 for(i=0;i<ns*ncatg;++i) p_lk[site*ncatgns + i] = 0.0;
1581 }
1582 }
1583 }
1584
1585 //////////////////////////////////////////////////////////////
1586 //////////////////////////////////////////////////////////////
1587
Default_Update_Partial_Lk(t_tree * tree,t_edge * b,t_node * d)1588 void Default_Update_Partial_Lk(t_tree *tree, t_edge *b, t_node *d)
1589 {
1590 /*
1591 |
1592 |<- b
1593 |
1594 d
1595 / \
1596 / \
1597 / \
1598 n_v1 n_v2
1599 */
1600 t_node *n_v1, *n_v2;//d's "left" and "right" neighbor nodes
1601 phydbl *p_lk;
1602 phydbl *p_lk_v1,*p_lk_v2;//Partial likelihood vector of node d, d's "left" neighbor, d's "right" neighbor. We fill *p_lk, and assume *p_lk_v1 and *p_lk_v2 are already filled.
1603 phydbl *Pij1,*Pij2;
1604 phydbl *tPij1,*tPij2;
1605 int *sum_scale, *sum_scale_v1, *sum_scale_v2;
1606 int *p_lk_loc;//Suppose site j, of a certain subtree, has "A" on one tip, and "C" on the other. If you come across this pattern again at site i<j, then you can simply copy the partial likelihoods
1607
1608 const unsigned int ncatg = tree->mod->ras->n_catg;
1609 const unsigned int ns = tree->mod->ns;
1610 const unsigned int n_patterns = tree->n_pattern;
1611
1612
1613 if(tree->n_root && tree->ignore_root == YES &&
1614 (d == tree->n_root->v[1] || d == tree->n_root->v[2]) &&
1615 (b == tree->n_root->b[1] || b == tree->n_root->b[2]))
1616 {
1617 assert(FALSE);
1618 }
1619
1620 if(d->tax)
1621 {
1622 PhyML_Fprintf(stderr,"\n. t_node %d is a leaf...",d->num);
1623 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s')\n",__FILE__,__LINE__,__FUNCTION__);
1624 Exit("\n");
1625 }
1626
1627
1628 n_v1 = n_v2 = NULL;
1629 p_lk = p_lk_v1 = p_lk_v2 = NULL;
1630 Pij1 = Pij2 = NULL;
1631 tPij1 = tPij2 = NULL;
1632 p_lk_loc = NULL;
1633 sum_scale_v1 = NULL;
1634 sum_scale_v2 = NULL;
1635
1636 Set_All_Partial_Lk(&n_v1,&n_v2,
1637 &p_lk,&sum_scale,&p_lk_loc,
1638 &Pij1,&tPij1,&p_lk_v1,&sum_scale_v1,
1639 &Pij2,&tPij2,&p_lk_v2,&sum_scale_v2,
1640 d,b,tree);
1641
1642 Core_Default_Update_Partial_Lk(n_v1,n_v2,
1643 p_lk,p_lk_v1,p_lk_v2,
1644 Pij1,Pij2,
1645 sum_scale,sum_scale_v1,sum_scale_v2,
1646 ns,ncatg,n_patterns,
1647 tree->apply_lk_scaling,
1648 tree->data->wght);
1649 }
1650
1651
1652 //////////////////////////////////////////////////////////////
1653 //////////////////////////////////////////////////////////////
1654
Core_Default_Update_Partial_Lk(const t_node * n_v1,const t_node * n_v2,phydbl * plk0,const phydbl * plk1,const phydbl * plk2,const phydbl * Pij1,const phydbl * Pij2,int * sum_scale0,const int * sum_scale1,const int * sum_scale2,const int ns,const int ncatg,const int npatterns,const int apply_scaling,const phydbl * wght)1655 void Core_Default_Update_Partial_Lk(const t_node *n_v1, const t_node *n_v2,
1656 phydbl *plk0, const phydbl *plk1, const phydbl *plk2,
1657 const phydbl *Pij1, const phydbl *Pij2,
1658 int *sum_scale0, const int *sum_scale1, const int *sum_scale2,
1659 const int ns, const int ncatg, const int npatterns, const int apply_scaling,
1660 const phydbl *wght)
1661 {
1662
1663 unsigned int i,site,ncatgns,catg,nsns;
1664 int state_v1,state_v2;
1665 int ambiguity_check_v1,ambiguity_check_v2;
1666 int sum_scale_v1_val, sum_scale_v2_val;
1667 phydbl largest_p_lk;
1668 const phydbl *init_Pij1, *init_Pij2;
1669
1670 ncatgns = ncatg*ns;
1671 nsns = ns*ns;
1672 init_Pij1 = Pij1;
1673 init_Pij2 = Pij2;
1674
1675 /* For every site in the alignment */
1676 for(site=0;site<npatterns;site++)
1677 {
1678 if(wght[site] > SMALL)
1679 {
1680 state_v1 = state_v2 = -1;
1681 ambiguity_check_v1 = ambiguity_check_v2 = YES;
1682
1683 /* n_v1 and n_v2 are tip nodes */
1684 if(n_v1->tax)
1685 {
1686 /* Is the state at this tip ambiguous? */
1687 ambiguity_check_v1 = n_v1->c_seq->is_ambigu[site];
1688 if(ambiguity_check_v1 == NO) state_v1 = n_v1->c_seq->d_state[site];
1689 }
1690
1691 if(n_v2->tax)
1692 {
1693 /* Is the state at this tip ambiguous? */
1694 ambiguity_check_v2 = n_v2->c_seq->is_ambigu[site];
1695 if(ambiguity_check_v2 == NO) state_v2 = n_v2->c_seq->d_state[site];
1696 }
1697
1698 Pij1 = init_Pij1;
1699 Pij2 = init_Pij2;
1700
1701 /* For all the rate classes */
1702 for(catg=0;catg<ncatg;++catg)
1703 {
1704 if(ambiguity_check_v1 == NO && ambiguity_check_v2 == NO)
1705 {
1706 Partial_Lk_Exex(Pij1,state_v1,
1707 Pij2,state_v2,
1708 ns,plk0);
1709 }
1710 else if(ambiguity_check_v1 == YES && ambiguity_check_v2 == NO)
1711 {
1712 Partial_Lk_Exin(Pij2,state_v2,
1713 Pij1,plk1,
1714 ns,plk0);
1715 }
1716 else if(ambiguity_check_v1 == NO && ambiguity_check_v2 == YES)
1717 {
1718 Partial_Lk_Exin(Pij1,state_v1,
1719 Pij2,plk2,
1720 ns,plk0);
1721 }
1722 else
1723 {
1724 Partial_Lk_Inin(Pij1,plk1,
1725 Pij2,plk2,
1726 ns,plk0);
1727 }
1728
1729 Pij1 += nsns;
1730 Pij2 += nsns;
1731
1732 plk1 += (n_v1->tax) ? 0 : ns;
1733 plk2 += (n_v2->tax) ? 0 : ns;
1734 plk0 += ns;
1735 }
1736
1737 plk1 += (n_v1->tax) ? ns : 0;
1738 plk2 += (n_v2->tax) ? ns : 0;
1739
1740 sum_scale_v1_val = (sum_scale1)?(sum_scale1[site]):(0);
1741 sum_scale_v2_val = (sum_scale2)?(sum_scale2[site]):(0);
1742 sum_scale0[site] = sum_scale_v1_val + sum_scale_v2_val;
1743
1744 plk0 -= ncatgns;
1745 largest_p_lk = -BIG;
1746 for(i=0;i<ncatgns;++i)
1747 if(plk0[i] > largest_p_lk)
1748 largest_p_lk = plk0[i];
1749
1750 if(largest_p_lk < INV_TWO_TO_THE_LARGE && apply_scaling == YES)
1751 {
1752 for(i=0;i<ncatgns;++i) plk0[i] *= TWO_TO_THE_LARGE;
1753 sum_scale0[site] += LARGE;
1754 }
1755 plk0 += ncatgns;
1756 }
1757 else
1758 {
1759 plk0 += ncatgns;
1760 plk1 += (n_v1->tax) ? ns : ncatgns;
1761 plk2 += (n_v2->tax) ? ns : ncatgns;
1762 }
1763 }
1764 }
1765 #endif
1766
1767 //////////////////////////////////////////////////////////////
1768 //////////////////////////////////////////////////////////////
1769
Return_Abs_Lk(t_tree * tree)1770 phydbl Return_Abs_Lk(t_tree *tree)
1771 {
1772 Lk(NULL,tree);
1773 return FABS(tree->c_lnL);
1774 }
1775
1776 //////////////////////////////////////////////////////////////
1777 //////////////////////////////////////////////////////////////
1778
ML_Dist(calign * data,t_mod * mod)1779 matrix *ML_Dist(calign *data, t_mod *mod)
1780 {
1781 int i,j,k,l;
1782 phydbl init;
1783 int n_catg;
1784 phydbl d_max,sum;
1785 matrix *mat;
1786 calign *twodata,*tmpdata;
1787 int state0, state1;
1788 phydbl *F,len;
1789 eigen *eigen_struct;
1790
1791 tmpdata = (calign *)mCalloc(1,sizeof(calign));
1792 tmpdata->c_seq = (align **)mCalloc(2,sizeof(align *));
1793 tmpdata->obs_state_frq = (phydbl *)mCalloc(mod->ns,sizeof(phydbl));
1794 tmpdata->ambigu = (short int *)mCalloc(data->crunch_len,sizeof(short int));
1795 F = (phydbl *)mCalloc(mod->ns*mod->ns,sizeof(phydbl ));
1796 eigen_struct = (eigen *)Make_Eigen_Struct(mod->ns);
1797
1798 Set_Update_Eigen(YES,mod);
1799 Update_Boundaries(mod);
1800 Update_Eigen(mod);
1801
1802 tmpdata->n_otu = 2;
1803
1804 tmpdata->crunch_len = data->crunch_len;
1805 tmpdata->init_len = data->init_len;
1806
1807 mat = NULL;
1808 if(mod->io->datatype == NT) mat = (mod->whichmodel < 10)?(K80_dist(data,1E+6)):(JC69_Dist(data,mod));
1809 else if(mod->io->datatype == AA) mat = JC69_Dist(data,mod);
1810 else if(mod->io->datatype == GENERIC) mat = JC69_Dist(data,mod);
1811
1812
1813 for(i=0;i<mod->ras->n_catg;i++) /* Don't use the discrete gamma distribution */
1814 {
1815 mod->ras->gamma_rr->v[i] = 1.0;
1816 mod->ras->gamma_r_proba->v[i] = 1.0;
1817 }
1818
1819 n_catg = mod->ras->n_catg;
1820 mod->ras->n_catg = 1;
1821
1822 for(j=0;j<data->n_otu-1;j++)
1823 {
1824 tmpdata->c_seq[0] = data->c_seq[j];
1825 tmpdata->c_seq[0]->name = data->c_seq[j]->name;
1826 tmpdata->wght = data->wght;
1827
1828 for(k=j+1;k<data->n_otu;k++)
1829 {
1830 tmpdata->c_seq[1] = data->c_seq[k];
1831 tmpdata->c_seq[1]->name = data->c_seq[k]->name;
1832
1833 twodata = Compact_Cdata(tmpdata,mod->io);
1834
1835 Check_Ambiguities(twodata,mod->io->datatype,mod->io->state_len);
1836
1837 Hide_Ambiguities(twodata);
1838
1839 init = mat->dist[j][k];
1840
1841 if((init > DIST_MAX-SMALL) || (init < .0)) init = 0.1;
1842
1843 d_max = init;
1844
1845 for(i=0;i<mod->ns*mod->ns;++i) F[i]=.0;
1846 len = 0.0;
1847 for(l=0;l<twodata->c_seq[0]->len;++l)
1848 {
1849 state0 = Assign_State(twodata->c_seq[0]->state+l*mod->io->state_len,mod->io->datatype,mod->io->state_len);
1850 state1 = Assign_State(twodata->c_seq[1]->state+l*mod->io->state_len,mod->io->datatype,mod->io->state_len);
1851
1852 if((state0 > -1) && (state1 > -1))
1853 {
1854 F[mod->ns*state0+state1] += twodata->wght[l];
1855 len += twodata->wght[l];
1856 }
1857 }
1858
1859 if(len > .0)
1860 {
1861 for(i=0;i<mod->ns*mod->ns;++i) F[i] /= len;
1862 }
1863
1864 sum = 0.;
1865 for(i=0;i<mod->ns*mod->ns;++i) sum += F[i];
1866
1867 /* for(i=0;i<mod->ns*mod->ns;++i) PhyML_Printf("\n. %g",F[i]); */
1868
1869 /* if(sum < .001) d_max = -1.; */
1870 if(sum < .001) d_max = init;
1871 else if((sum > 1. - .001) && (sum < 1. + .001)) Opt_Dist_F(&(d_max),F,mod);
1872 else
1873 {
1874 PhyML_Fprintf(stderr,"\n\n. sum = %f",sum);
1875 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1876 Exit("");
1877 }
1878
1879 if(d_max >= DIST_MAX) d_max = DIST_MAX;
1880
1881
1882 /* Do not correct for dist < BL_MIN, otherwise Fill_Missing_Dist
1883 * will not be called
1884 */
1885 mat->dist[j][k] = d_max;
1886 mat->dist[k][j] = mat->dist[j][k];
1887 Free_Calign(twodata);
1888 }
1889 }
1890
1891 mod->ras->n_catg = n_catg;
1892
1893
1894 Free(tmpdata->ambigu);
1895 Free(tmpdata->obs_state_frq);
1896 Free(tmpdata->c_seq);
1897 free(tmpdata);
1898 Free_Eigen(eigen_struct);
1899 Free(F);
1900
1901 return mat;
1902 }
1903
1904 //////////////////////////////////////////////////////////////
1905 //////////////////////////////////////////////////////////////
1906
Lk_Given_Two_Seq(calign * data,int numseq1,int numseq2,phydbl dist,t_mod * mod,phydbl * loglk)1907 phydbl Lk_Given_Two_Seq(calign *data, int numseq1, int numseq2, phydbl dist, t_mod *mod, phydbl *loglk)
1908 {
1909 align *seq1,*seq2;
1910 phydbl site_lk,log_site_lk;
1911 int i,j,k,l;
1912 /* phydbl **p_lk_l,**p_lk_r; */
1913 phydbl *p_lk_l,*p_lk_r;
1914 phydbl len;
1915 int dim1,dim2;
1916
1917 dim1 = mod->ns;
1918 dim2 = mod->ns * mod->ns;
1919
1920 DiscreteGamma(mod->ras->gamma_r_proba->v, mod->ras->gamma_rr->v, mod->ras->alpha->v,
1921 mod->ras->alpha->v,mod->ras->n_catg,mod->ras->gamma_median);
1922
1923 seq1 = data->c_seq[numseq1];
1924 seq2 = data->c_seq[numseq2];
1925
1926
1927 p_lk_l = (phydbl *)mCalloc(data->c_seq[0]->len * mod->ns,sizeof(phydbl));
1928 p_lk_r = (phydbl *)mCalloc(data->c_seq[0]->len * mod->ns,sizeof(phydbl));
1929
1930
1931 for(i=0;i<mod->ras->n_catg;i++)
1932 {
1933 len = dist*mod->ras->gamma_rr->v[i];
1934 if(len < mod->l_min) len = mod->l_min;
1935 else if(len > mod->l_max) len = mod->l_max;
1936 PMat(len,mod,dim2*i,mod->Pij_rr->v,NULL);
1937 }
1938
1939 if(mod->io->datatype == NT)
1940 {
1941 For(i,data->c_seq[0]->len)
1942 {
1943 Init_Tips_At_One_Site_Nucleotides_Float(seq1->state[i],i*mod->ns,p_lk_l);
1944 Init_Tips_At_One_Site_Nucleotides_Float(seq2->state[i],i*mod->ns,p_lk_r);
1945 }
1946 }
1947 else if(mod->io->datatype == AA)
1948 {
1949 For(i,data->c_seq[0]->len)
1950 {
1951 Init_Tips_At_One_Site_AA_Float(seq1->state[i],i*mod->ns,p_lk_l);
1952 Init_Tips_At_One_Site_AA_Float(seq2->state[i],i*mod->ns,p_lk_r);
1953 }
1954 }
1955 else
1956 {
1957 PhyML_Fprintf(stderr,"\n\n. Not implemented yet...");
1958 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1959 Warn_And_Exit("\n");
1960 }
1961
1962
1963 site_lk = .0;
1964 *loglk = 0;
1965
1966 For(i,data->c_seq[0]->len)
1967 {
1968 if(data->wght[i] > 0.0)
1969 {
1970 site_lk = log_site_lk = .0;
1971 if(!data->ambigu[i])
1972 {
1973 for(k=0;k<mod->ns;k++) {if(p_lk_l[i*mod->ns+k] > .0001) break;}
1974 for(l=0;l<mod->ns;l++) {if(p_lk_r[i*mod->ns+l] > .0001) break;}
1975 for(j=0;j<mod->ras->n_catg;j++)
1976 {
1977 site_lk +=
1978 mod->ras->gamma_r_proba->v[j] *
1979 mod->e_frq->pi->v[k] *
1980 p_lk_l[i*dim1+k] *
1981 mod->Pij_rr->v[j*dim2+k*dim1+l] *
1982 p_lk_r[i*dim1+l];
1983 }
1984 }
1985 else
1986 {
1987 for(j=0;j<mod->ras->n_catg;j++)
1988 {
1989 for(k=0;k<mod->ns;k++) /*sort sum terms ? No global effect*/
1990 {
1991 for(l=0;l<mod->ns;l++)
1992 {
1993 site_lk +=
1994 mod->ras->gamma_r_proba->v[j] *
1995 mod->e_frq->pi->v[k] *
1996 p_lk_l[i*dim1+k] *
1997 mod->Pij_rr->v[j*dim2+k*dim1+l] *
1998 p_lk_r[i*dim1+l];
1999 }
2000 }
2001 }
2002 }
2003
2004 if(site_lk <= .0)
2005 {
2006 PhyML_Fprintf(stderr,"\n\n. '%c' '%c'\n",seq1->state[i],seq2->state[i]);
2007 Exit("\n. Err: site lk <= 0\n");
2008 }
2009
2010 log_site_lk += (phydbl)log(site_lk);
2011
2012 *loglk += data->wght[i] * log_site_lk;/* sort sum terms ? No global effect*/
2013 }
2014 }
2015
2016 /* For(i,data->c_seq[0]->len) */
2017 /* { */
2018 /* Free(p_lk_l[i]); */
2019 /* Free(p_lk_r[i]); */
2020 /* } */
2021
2022 Free(p_lk_l); Free(p_lk_r);
2023 return *loglk;
2024 }
2025
2026 //////////////////////////////////////////////////////////////
2027 //////////////////////////////////////////////////////////////
2028 // Multinomial log likelihood
Unconstraint_Lk(t_tree * tree)2029 void Unconstraint_Lk(t_tree *tree)
2030 {
2031 int i;
2032
2033 tree->unconstraint_lk = .0;
2034 for(i=0;i<tree->data->crunch_len;i++) tree->unconstraint_lk += tree->data->wght[i]*(phydbl)log(tree->data->wght[i]);
2035 tree->unconstraint_lk -= tree->data->init_len*(phydbl)log(tree->data->init_len);
2036 }
2037
2038 //////////////////////////////////////////////////////////////
2039 //////////////////////////////////////////////////////////////
2040 // Log-likelihood assuming a tree with edge of infinite lengths
Composite_Lk(t_tree * tree)2041 void Composite_Lk(t_tree *tree)
2042 {
2043 int i;
2044 tree->composite_lk = 0.0;
2045 for(i=0;i<tree->mod->ns;++i)
2046 tree->composite_lk +=
2047 tree->data->obs_state_frq[i]*
2048 tree->data->init_len*
2049 tree->n_otu*
2050 log(tree->mod->e_frq->pi->v[i]);
2051 }
2052
2053 //////////////////////////////////////////////////////////////
2054 //////////////////////////////////////////////////////////////
2055
Init_Partial_Lk_Tips_Double(t_tree * tree)2056 void Init_Partial_Lk_Tips_Double(t_tree *tree)
2057 {
2058 unsigned int curr_site,i,dim1;
2059
2060 if(tree->is_mixt_tree == YES) return;
2061
2062 dim1 = tree->mod->ns;
2063
2064
2065 for(i=0;i<tree->n_otu;i++)
2066 {
2067 if(!tree->a_nodes[i]->c_seq ||
2068 strcmp(tree->a_nodes[i]->c_seq->name,tree->a_nodes[i]->name))
2069 {
2070 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
2071 Exit("");
2072 }
2073 }
2074
2075 for(curr_site=0;curr_site<tree->data->crunch_len;curr_site++)
2076 {
2077 for(i=0;i<tree->n_otu;i++)
2078 {
2079 if (tree->io->datatype == NT)
2080 Init_Tips_At_One_Site_Nucleotides_Float(tree->a_nodes[i]->c_seq->state[curr_site],
2081 curr_site*dim1,
2082 tree->a_nodes[i]->b[0]->p_lk_tip_r);
2083 else if(tree->io->datatype == AA)
2084 Init_Tips_At_One_Site_AA_Float(tree->a_nodes[i]->c_seq->state[curr_site],
2085 curr_site*dim1,
2086 tree->a_nodes[i]->b[0]->p_lk_tip_r);
2087
2088 else if(tree->io->datatype == GENERIC)
2089 Init_Tips_At_One_Site_Generic_Float(tree->a_nodes[i]->c_seq->state+curr_site*tree->mod->io->state_len,
2090 tree->mod->ns,
2091 tree->mod->io->state_len,
2092 curr_site*dim1,
2093 tree->a_nodes[i]->b[0]->p_lk_tip_r);
2094 }
2095 }
2096 }
2097
2098 //////////////////////////////////////////////////////////////
2099 //////////////////////////////////////////////////////////////
2100
Init_Partial_Lk_Tips_Int(t_tree * tree)2101 void Init_Partial_Lk_Tips_Int(t_tree *tree)
2102 {
2103 int curr_site,i,dim1;
2104
2105 if(tree->is_mixt_tree == YES) return;
2106
2107 dim1 = tree->mod->ns;
2108
2109 for(i=0;i<tree->n_otu;i++)
2110 {
2111 if(!tree->a_nodes[i]->c_seq ||
2112 strcmp(tree->a_nodes[i]->c_seq->name,tree->a_nodes[i]->name))
2113 {
2114 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
2115 Exit("");
2116 }
2117 }
2118
2119 for(curr_site=0;curr_site<tree->data->crunch_len;curr_site++)
2120 {
2121 for(i=0;i<tree->n_otu;i++)
2122 {
2123 /* printf("\n. site: %3d %c",curr_site,tree->a_nodes[i]->c_seq->state[curr_site]); */
2124 /* printf("\n. init at %s %p",tree->a_nodes[i]->name,tree->a_nodes[i]->b[0]->p_lk_tip_r); fflush(NULL); */
2125 if(tree->io->datatype == NT)
2126 {
2127 Init_Tips_At_One_Site_Nucleotides_Float(tree->a_nodes[i]->c_seq->state[curr_site],
2128 curr_site*dim1,
2129 tree->a_nodes[i]->b[0]->p_lk_tip_r);
2130 /* Init_Tips_At_One_Site_Nucleotides_Int(tree->data->c_seq[i]->state[curr_site], */
2131 /* curr_site*dim1, */
2132 /* tree->a_nodes[i]->b[0]->p_lk_tip_r); */
2133 }
2134 else if(tree->io->datatype == AA)
2135 Init_Tips_At_One_Site_AA_Float(tree->a_nodes[i]->c_seq->state[curr_site],
2136 curr_site*dim1,
2137 tree->a_nodes[i]->b[0]->p_lk_tip_r);
2138 /* Init_Tips_At_One_Site_AA_Int(tree->data->c_seq[i]->state[curr_site], */
2139 /* curr_site*dim1, */
2140 /* tree->a_nodes[i]->b[0]->p_lk_tip_r); */
2141
2142 else if(tree->io->datatype == GENERIC)
2143 {
2144 Init_Tips_At_One_Site_Generic_Float(tree->a_nodes[i]->c_seq->state+curr_site*tree->mod->io->state_len,
2145 tree->mod->ns,
2146 tree->mod->io->state_len,
2147 curr_site*dim1,
2148 tree->a_nodes[i]->b[0]->p_lk_tip_r);
2149
2150 /* Init_Tips_At_One_Site_Generic_Int(tree->data->c_seq[i]->state+curr_site*tree->mod->io->state_len, */
2151 /* tree->mod->ns, */
2152 /* tree->mod->io->state_len, */
2153 /* curr_site*dim1, */
2154 /* tree->a_nodes[i]->b[0]->p_lk_tip_r); */
2155 }
2156 #ifdef BEAGLE
2157 //Recall that tip partials are stored on the branch leading
2158 //to the tip, rather than on the tip itself (hence `p_lk_tip_idx`
2159 //is a field of the branch (i.e. b[0]) rather than the node.
2160 //Secondly, the BEAGLE's partial buffers are laid out as
2161 //BEAGLE's partials buffer = [ tax1, tax2, ..., taxN, b1Left, b2Left, b3Left,...,bMLeft, b1Rght, b2Rght, b3Rght,...,bMRght] (N taxa, M branches)
2162 tree->a_nodes[i]->b[0]->p_lk_tip_idx = i;
2163 #endif
2164 }
2165 }
2166 }
2167
2168 //////////////////////////////////////////////////////////////
2169 //////////////////////////////////////////////////////////////
2170
Update_PMat_At_Given_Edge(t_edge * b_fcus,t_tree * tree)2171 void Update_PMat_At_Given_Edge(t_edge *b_fcus, t_tree *tree)
2172 {
2173 int i;
2174 phydbl len;
2175 phydbl l_min, l_max;
2176 phydbl shape, scale, mean, var;
2177
2178 assert(b_fcus);
2179 assert(tree);
2180
2181 if(tree->is_mixt_tree == YES)
2182 {
2183 MIXT_Update_PMat_At_Given_Edge(b_fcus,tree);
2184 return;
2185 }
2186
2187 if(tree->mixt_tree != NULL) assert(tree->mod->ras->n_catg == 1);
2188
2189 if(tree->io->mod->gamma_mgf_bl == YES) Set_Br_Len_Var(b_fcus,tree);
2190
2191 l_min = tree->mod->l_min;
2192 l_max = tree->mod->l_max;
2193
2194
2195 len = -1.0;
2196
2197 if(tree->mod->log_l == YES) b_fcus->l->v = exp(b_fcus->l->v);
2198
2199 /* if(b_fcus->l->v < l_min) b_fcus->l->v = l_min; */
2200 /* if(b_fcus->l->v > l_max) b_fcus->l->v = l_max; */
2201
2202 for(i=0;i<tree->mod->ras->n_catg;i++)
2203 {
2204 if(tree->mod->ras->skip_rate_cat[i] == YES) continue;
2205
2206 //Update the branch length
2207 if(b_fcus->has_zero_br_len == YES)
2208 {
2209 #ifdef BEAGLE
2210 Warn_And_Exit(TODO_BEAGLE);
2211 #endif
2212 len = -1.0;
2213 mean = -1.0;
2214 var = -1.0;
2215 }
2216 else
2217 {
2218 len = MAX(0.0,b_fcus->l->v)*tree->mod->ras->gamma_rr->v[i];//branch_len * rate
2219 len *= tree->mod->br_len_mult->v;
2220 if(tree->mixt_tree) len *= tree->mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number];
2221 if(len < l_min) len = l_min;
2222 else if(len > l_max) len = l_max;
2223
2224 mean = len;
2225 var = MAX(0.0,b_fcus->l_var->v) * POW(tree->mod->ras->gamma_rr->v[i]*tree->mod->br_len_mult->v,2);
2226 if(tree->mixt_tree) var *= POW(tree->mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number],2);
2227
2228 /* var = 1.E-10; */
2229
2230 if(var > tree->mod->l_var_max) var = tree->mod->l_var_max;
2231 if(var < tree->mod->l_var_min) var = tree->mod->l_var_min;
2232 }
2233
2234 //Update the transition prob. matrix
2235 if(tree->mod->gamma_mgf_bl == NO)
2236 {
2237 #ifdef BEAGLE
2238 assert(UNINITIALIZED != tree->mod->b_inst);
2239 #endif
2240
2241 PMat(len,tree->mod,i*tree->mod->ns*tree->mod->ns,b_fcus->Pij_rr,b_fcus->tPij_rr);
2242 }
2243 else
2244 {
2245 #ifdef BEAGLE
2246 Warn_And_Exit(TODO_BEAGLE);
2247 #endif
2248
2249 shape = mean*mean/var;
2250 scale = var/mean;
2251 PMat_MGF_Gamma(b_fcus->Pij_rr+tree->mod->ns*tree->mod->ns*i,shape,scale,1.0,tree->mod);
2252 }
2253 }
2254
2255 #ifdef BEAGLE
2256 int whichmodel = tree->mod->whichmodel;
2257 //Only for some models we use Beagle to compute/update the P-matrices, for other models
2258 //we compute them in PhyML and explicitly set the P-matrices in BEAGLE
2259 if((tree->mod->io->datatype == AA || whichmodel==GTR || whichmodel==CUSTOM) && tree->mod->use_m4mod == NO)
2260 {
2261 if(b_fcus->has_zero_br_len == YES)
2262 Warn_And_Exit(TODO_BEAGLE);
2263
2264 //
2265 update_beagle_eigen(tree->mod);
2266 update_beagle_ras(tree->mod);
2267
2268 //
2269 len = MAX(0.0, b_fcus->l->v) * tree->mod->br_len_mult->v;
2270 int p_matrices[1] = b_fcus->Pij_rr_idx;
2271 double branch_lens[1] = len;
2272 int ret = beagleUpdateTransitionMatrices(tree->b_inst,0,p_matrices,NULL,NULL,branch_lens,1);
2273 if(ret<0)
2274 {
2275 PhyML_Fprintf(stderr, "beagleUpdateTransitionMatrices() on instance %i failed:%i\n\n",tree->b_inst,ret);
2276 Exit("");
2277 }
2278 //Retrieve a "local" copy of the P-matrix
2279 ret = beagleGetTransitionMatrix(tree->b_inst, b_fcus->Pij_rr_idx, b_fcus->Pij_rr);
2280 if(ret<0)
2281 {
2282 PhyML_Fprintf(stderr, "beagleGetTransitionMatrix() on instance %i failed:%i\n\n",tree->b_inst,ret);
2283 Exit("");
2284 }
2285 }
2286 else
2287 {
2288 int ret = beagleSetTransitionMatrix(tree->b_inst, b_fcus->Pij_rr_idx, b_fcus->Pij_rr, -1);
2289 if(ret<0)
2290 {
2291 PhyML_Fprintf(stderr, "beagleSetTransitionMatrix() on instance %i failed:%i\n\n",tree->b_inst,ret);
2292 Exit("");
2293 }
2294 }
2295 #endif
2296 if(tree->mod->log_l == YES) b_fcus->l->v = log(b_fcus->l->v);
2297
2298 // Print_Model(tree->mod);
2299 // Dump_Arr_D(tree->cur_site_lk, tree->n_pattern);
2300 // Print_Edge_PMats(tree, b_fcus);
2301 // Print_Edge_Likelihoods(tree,b_fcus,true);
2302 }
2303
2304 //////////////////////////////////////////////////////////////
2305 //////////////////////////////////////////////////////////////
2306
Update_Partial_Lk_Along_A_Path(t_node ** path,int path_length,t_tree * tree)2307 void Update_Partial_Lk_Along_A_Path(t_node **path, int path_length, t_tree *tree)
2308 {
2309 int i,j;
2310
2311 for(i=0;i<path_length-1;++i)
2312 {
2313 for(j=0;j<3;++j)
2314 if(path[i]->v[j] == path[i+1])
2315 {
2316 if(path[i] == path[i]->b[j]->left)
2317 {
2318 Update_Partial_Lk(tree,path[i]->b[j],path[i]->b[j]->left);
2319 }
2320 else if(path[i] == path[i]->b[j]->rght)
2321 {
2322 Update_Partial_Lk(tree,path[i]->b[j],path[i]->b[j]->rght);
2323 }
2324 else
2325 {
2326 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d. \n",__FILE__,__LINE__);
2327 assert(FALSE);
2328 }
2329 break;
2330 }
2331 #ifdef DEBUG
2332 if(j == 3)
2333 {
2334 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d.\n",__FILE__,__LINE__);
2335 assert(FALSE);
2336 }
2337 #endif
2338 }
2339 }
2340
2341 //////////////////////////////////////////////////////////////
2342 //////////////////////////////////////////////////////////////
2343
Lk_Dist(phydbl * F,phydbl dist,t_mod * mod)2344 phydbl Lk_Dist(phydbl *F, phydbl dist, t_mod *mod)
2345 {
2346 int i,j,k;
2347 phydbl lnL,len;
2348 int dim1,dim2;
2349 phydbl pi, pijk;
2350
2351 // Compute likelihood of the model made of the
2352 // first class of the mixture.
2353 /* if(mod->is_mixt_mod == YES) mod = mod->next; */
2354 /* assert(mod); */
2355
2356 if(mod->log_l == YES) dist = exp(dist);
2357
2358 for(k=0;k<mod->ras->n_catg;k++)
2359 {
2360 len = dist*mod->ras->gamma_rr->v[k];
2361 if(len < mod->l_min) len = mod->l_min;
2362 else if(len > mod->l_max) len = mod->l_max;
2363 PMat(len,mod,mod->ns*mod->ns*k,mod->Pij_rr->v,NULL);
2364 /* PhyML_Printf("\n. p: %g len: %g",mod->Pij_rr->v[0],len); */
2365 }
2366
2367 dim1 = mod->ns*mod->ns;
2368 dim2 = mod->ns;
2369 lnL = .0;
2370 pi = -1.;
2371 pijk = -1.;
2372
2373 for(i=0;i<mod->ns-1;i++)
2374 {
2375 pi = mod->e_frq->pi->v[i];
2376
2377 for(j=i+1;j<mod->ns;j++)
2378 {
2379 for(k=0;k<mod->ras->n_catg;k++)
2380 {
2381 pijk = mod->Pij_rr->v[dim1*k+dim2*i+j];
2382 lnL += (F[dim1*k+dim2*i+j] + F[dim1*k+dim2*j+i]) * log(pi * pijk);
2383 /* PhyML_Printf("\n. pijk: %g lnL: %g",pijk,lnL); */
2384 }
2385 }
2386 }
2387
2388 for(i=0;i<mod->ns;i++)
2389 {
2390 pi = mod->e_frq->pi->v[i];
2391
2392 for(k=0;k<mod->ras->n_catg;k++)
2393 {
2394 pijk = mod->Pij_rr->v[dim1*k+dim2*i+i];
2395 lnL += F[dim1*k+dim2*i+i]* log(pi * pijk);
2396 /* PhyML_Printf("\n. pijk: %g lnL: %g",pijk,lnL); */
2397 }
2398 }
2399
2400 return lnL;
2401 }
2402
2403 //////////////////////////////////////////////////////////////
2404 //////////////////////////////////////////////////////////////
2405
Update_Lk_At_Given_Edge(t_edge * b_fcus,t_tree * tree)2406 phydbl Update_Lk_At_Given_Edge(t_edge *b_fcus, t_tree *tree)
2407 {
2408 Update_Partial_Lk(tree,b_fcus,b_fcus->left);
2409 Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
2410 tree->c_lnL = Lk(b_fcus,tree);
2411 return tree->c_lnL;
2412 }
2413
2414 //////////////////////////////////////////////////////////////
2415 //////////////////////////////////////////////////////////////
2416
Init_Partial_Lk_Loc(t_tree * tree)2417 void Init_Partial_Lk_Loc(t_tree *tree)
2418 {
2419 int i,j;
2420 t_node *d;
2421 int *patt_id_d;
2422
2423 if(tree->is_mixt_tree == YES) return;
2424
2425 for(i=0;i<2*tree->n_otu-1;++i)
2426 {
2427 for(j=0;j<tree->n_pattern;j++)
2428 {
2429 tree->a_edges[i]->p_lk_loc_left[j] = j;
2430 tree->a_edges[i]->p_lk_loc_rght[j] = j;
2431 }
2432 }
2433
2434 for(i=0;i<tree->n_otu;i++)
2435 {
2436 d = tree->a_nodes[i];
2437 patt_id_d = (d == d->b[0]->left)?(d->b[0]->patt_id_left):(d->b[0]->patt_id_rght);
2438 for(j=0;j<tree->n_pattern;j++)
2439 {
2440 assert(tree->a_nodes[d->num]->c_seq);
2441 patt_id_d[j] = (int)tree->a_nodes[d->num]->c_seq->state[j];
2442 }
2443 }
2444 }
2445
2446 //////////////////////////////////////////////////////////////
2447 //////////////////////////////////////////////////////////////
2448
Lk_Normal_Approx(t_tree * tree)2449 phydbl Lk_Normal_Approx(t_tree *tree)
2450 {
2451 phydbl lnL;
2452 int i;
2453 int dim;
2454 phydbl first_order;
2455
2456 dim = 2*tree->n_otu-3;
2457
2458 lnL = Dnorm_Multi_Given_InvCov_Det(tree->rates->u_cur_l,
2459 tree->rates->mean_l,
2460 tree->rates->invcov,
2461 tree->rates->covdet,
2462 2*tree->n_otu-3,YES);
2463
2464 first_order = 0.0;
2465 for(i=0;i<dim;i++) first_order += (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]) * tree->rates->grad_l[i];
2466
2467 lnL += first_order;
2468
2469 return(lnL);
2470
2471 }
2472
2473 //////////////////////////////////////////////////////////////
2474 //////////////////////////////////////////////////////////////
2475
2476
Wrap_Part_Lk_At_Given_Edge(t_edge * b,t_tree * tree,supert_tree * stree)2477 phydbl Wrap_Part_Lk_At_Given_Edge(t_edge *b, t_tree *tree, supert_tree *stree)
2478 {
2479 return -1.0;
2480 /* return PART_Lk_At_Given_Edge(b,stree);; */
2481 }
2482
2483 //////////////////////////////////////////////////////////////
2484 //////////////////////////////////////////////////////////////
2485
Wrap_Part_Lk(t_edge * b,t_tree * tree,supert_tree * stree)2486 phydbl Wrap_Part_Lk(t_edge *b, t_tree *tree, supert_tree *stree)
2487 {
2488 return -1.0;
2489 /* return PART_Lk(stree); */
2490 }
2491
2492 //////////////////////////////////////////////////////////////
2493 //////////////////////////////////////////////////////////////
2494
2495
Wrap_Lk(t_edge * b,t_tree * tree,supert_tree * stree)2496 phydbl Wrap_Lk(t_edge *b, t_tree *tree, supert_tree *stree)
2497 {
2498 Lk(NULL,tree);
2499 return tree->c_lnL;
2500 }
2501
2502 //////////////////////////////////////////////////////////////
2503 //////////////////////////////////////////////////////////////
2504
2505
Wrap_Geo_Lk(t_edge * b,t_tree * tree,supert_tree * stree)2506 phydbl Wrap_Geo_Lk(t_edge *b, t_tree *tree, supert_tree *stree)
2507 {
2508 TIPO_Lk(tree);
2509 return tree->geo_lnL;
2510 }
2511
2512 //////////////////////////////////////////////////////////////
2513 //////////////////////////////////////////////////////////////
2514
2515
Wrap_Lk_At_Given_Edge(t_edge * b,t_tree * tree,supert_tree * stree)2516 phydbl Wrap_Lk_At_Given_Edge(t_edge *b, t_tree *tree, supert_tree *stree)
2517 {
2518 Lk(b,tree);
2519 return tree->c_lnL;
2520 }
2521
2522 //////////////////////////////////////////////////////////////
2523 //////////////////////////////////////////////////////////////
2524
2525
Wrap_Lk_Rates(t_edge * b,t_tree * tree,supert_tree * stree)2526 phydbl Wrap_Lk_Rates(t_edge *b, t_tree *tree, supert_tree *stree)
2527 {
2528 RATES_Lk_Rates(tree);
2529 return tree->rates->c_lnL_rates;
2530 }
2531
2532 //////////////////////////////////////////////////////////////
2533 //////////////////////////////////////////////////////////////
2534
Wrap_Lk_Times(t_edge * b,t_tree * tree,supert_tree * stree)2535 phydbl Wrap_Lk_Times(t_edge *b, t_tree *tree, supert_tree *stree)
2536 {
2537 TIMES_Lk_Times(NO,tree);
2538 return tree->times->c_lnL_times;
2539 }
2540
2541 //////////////////////////////////////////////////////////////
2542 //////////////////////////////////////////////////////////////
2543
Wrap_Lk_Linreg(t_edge * b,t_tree * tree,supert_tree * stree)2544 phydbl Wrap_Lk_Linreg(t_edge *b, t_tree *tree, supert_tree *stree)
2545 {
2546 /* RATES_Lk_Linreg(tree); */
2547 return -1.;
2548 /* return tree->rates->c_lnL_linreg; */
2549 }
2550
2551 //////////////////////////////////////////////////////////////
2552 //////////////////////////////////////////////////////////////
2553
2554
Wrap_Diff_Lk_Norm_At_Given_Edge(t_edge * b,t_tree * tree,supert_tree * stree)2555 phydbl Wrap_Diff_Lk_Norm_At_Given_Edge(t_edge *b, t_tree *tree, supert_tree *stree)
2556 {
2557 phydbl diff;
2558 diff = Diff_Lk_Norm_At_Given_Edge(b,tree);
2559 return(-diff);
2560
2561 }
2562
2563 //////////////////////////////////////////////////////////////
2564 //////////////////////////////////////////////////////////////
2565
Check_Lk_At_Given_Edge(int verbose,t_tree * tree)2566 int Check_Lk_At_Given_Edge(int verbose, t_tree *tree)
2567 {
2568 int res;
2569 int i;
2570 phydbl *lk;
2571
2572 lk = (phydbl *)mCalloc(2*tree->n_otu-3,sizeof(phydbl));
2573
2574 res = 0;
2575 for(i=0;i<2*tree->n_otu-3;++i)
2576 {
2577 lk[i] = Lk(tree->a_edges[i],tree);
2578 if(verbose == YES) PhyML_Printf("\n. Edge %3d %13G %f %13G",
2579 tree->a_edges[i]->num,tree->a_edges[i]->l->v,lk[i],
2580 tree->a_edges[i]->l_var->v);
2581
2582 }
2583
2584 if(tree->n_root && tree->ignore_root == NO)
2585 {
2586 Lk(tree->n_root->b[1],tree);
2587 if(verbose == YES) PhyML_Printf("\nx Edge %3d %13G %f %13G",
2588 tree->n_root->b[1]->num,tree->n_root->b[1]->l->v,tree->c_lnL,
2589 tree->n_root->b[1]->l_var->v
2590 );
2591
2592 Lk(tree->n_root->b[2],tree);
2593 if(verbose == YES) PhyML_Printf("\nx Edge %3d %13G %f %13G",
2594 tree->n_root->b[2]->num,tree->n_root->b[2]->l->v,tree->c_lnL,
2595 tree->n_root->b[2]->l_var->v
2596 );
2597
2598 }
2599
2600 res=1;
2601 for(i=1;i<2*tree->n_otu-3;i++)
2602 {
2603 if(FABS(lk[i]-lk[i-1]) > 1.E-2) res=0;
2604 }
2605 Free(lk);
2606
2607 return res;
2608 }
2609
2610 //////////////////////////////////////////////////////////////
2611 //////////////////////////////////////////////////////////////
2612
2613 //////////////////////////////////////////////////////////////
2614 //////////////////////////////////////////////////////////////
2615 // Computes the value of fact_sum_scale, the part of the scaling factors
2616 // that is common to all classes of the mixture and scale the
2617 // likelihood for each mixture using the part of the scaling
2618 // factors that is class-specific
2619
Pull_Scaling_Factors(int site,t_edge * b,t_tree * tree)2620 void Pull_Scaling_Factors(int site, t_edge *b, t_tree *tree)
2621 {
2622 unsigned int catg;
2623 const unsigned int ncatg = tree->mod->ras->n_catg;
2624
2625 if(tree->apply_lk_scaling == NO)
2626 {
2627 tree->fact_sum_scale[site] = 0;
2628 for(catg=0;catg<ncatg;++catg) tree->unscaled_site_lk_cat[site*ncatg+catg] = tree->site_lk_cat[catg];
2629 return;
2630 }
2631 else
2632 {
2633 switch(tree->scaling_method)
2634 {
2635 case SCALE_RATE_SPECIFIC :
2636 {
2637 int *sum_scale_left_cat,*sum_scale_rght_cat;
2638 int exponent;
2639 phydbl max_sum_scale,min_sum_scale;
2640 phydbl sum,tmp,dum;
2641
2642 sum_scale_left_cat = b->sum_scale_left_cat;
2643 sum_scale_rght_cat = b->sum_scale_rght_cat;
2644
2645 max_sum_scale = (phydbl)BIG;
2646 min_sum_scale = -(phydbl)BIG;
2647
2648 for(catg=0;catg<ncatg;++catg)
2649 {
2650 sum_scale_left_cat[catg] =
2651 (b->sum_scale_left)?
2652 (b->sum_scale_left[site*ncatg+catg]):
2653 (0.0);
2654
2655 sum_scale_rght_cat[catg] =
2656 (b->sum_scale_rght)?
2657 (b->sum_scale_rght[site*ncatg+catg]):
2658 (0.0);
2659
2660 sum = sum_scale_left_cat[catg] + sum_scale_rght_cat[catg];
2661
2662 if(sum < .0)
2663 {
2664 PhyML_Fprintf(stderr,"\n. tree: %s\n",Write_Tree(tree));
2665 PhyML_Fprintf(stderr,"\n. b->num = %d sum = %G root ? %d",sum,b->num,b == tree->e_root);
2666 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d.\n",__FILE__,__LINE__);
2667 Exit("\n");
2668 }
2669
2670 dum = log(FABS(tree->site_lk_cat[catg]));
2671
2672 tmp = sum + ((phydbl)LOGBIG - dum)/(phydbl)LOG2;
2673 if(tmp < max_sum_scale) max_sum_scale = tmp; /* min of the maxs */
2674
2675 tmp = sum + ((phydbl)LOGSMALL - dum)/(phydbl)LOG2;
2676 if(tmp > min_sum_scale) min_sum_scale = tmp; /* max of the mins */
2677
2678 assert(isnan(tmp) == NO);
2679 }
2680
2681 if(min_sum_scale > max_sum_scale)
2682 {
2683 #ifdef SAFEMODE
2684 PhyML_Printf("\n. Numerical precision issue alert.");
2685 PhyML_Printf("\n. min_sum_scale = %G max_sum_scale = %G",min_sum_scale,max_sum_scale);
2686 #endif
2687 min_sum_scale = max_sum_scale;
2688 }
2689
2690 tree->fact_sum_scale[site] = (int)((max_sum_scale + min_sum_scale) / 2);
2691
2692 /* Apply scaling factors */
2693 for(catg=0;catg<ncatg;++catg)
2694 {
2695 exponent = -(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+tree->fact_sum_scale[site];
2696 Rate_Correction(exponent,tree->site_lk_cat + catg);
2697 }
2698
2699 break;
2700 }
2701 case SCALE_FAST :
2702 {
2703 int sum_scale_left,sum_scale_rght;
2704
2705 sum_scale_left =
2706 (b->sum_scale_left)?
2707 (b->sum_scale_left[site]):
2708 (0.0);
2709
2710 sum_scale_rght =
2711 (b->sum_scale_rght)?
2712 (b->sum_scale_rght[site]):
2713 (0.0);
2714
2715 tree->fact_sum_scale[site] = sum_scale_left + sum_scale_rght;
2716
2717 break;
2718 }
2719 }
2720 for(catg=0;catg<ncatg;++catg) tree->unscaled_site_lk_cat[site*ncatg+catg] = tree->site_lk_cat[catg];
2721 }
2722 }
2723
2724 //////////////////////////////////////////////////////////////
2725 //////////////////////////////////////////////////////////////
2726
2727 // Tree should be ready for likelihood analysis when calling
2728 // this function.
Stepwise_Add_Lk(t_tree * tree)2729 void Stepwise_Add_Lk(t_tree *tree)
2730 {
2731 t_edge **residuals,**targets,*best_target;
2732 int *nd_idx,i,j,n_targets,*tg_idx,n_opt;
2733
2734 residuals = (t_edge **)mCalloc(tree->n_otu-3,sizeof(t_edge *));
2735 targets = (t_edge **)mCalloc(2*tree->n_otu-3,sizeof(t_edge *));
2736 best_target = NULL;
2737 nd_idx = Permutate(tree->n_otu-3);
2738
2739 // Remove all tips except that corresponding to a_nodes[0],
2740 // a_nodes[1] and a_nodes[2].
2741 for(i=0;i<tree->n_otu-3;i++)
2742 {
2743 Prune_Subtree(tree->a_nodes[i+3]->v[0],
2744 tree->a_nodes[i+3],
2745 NULL,
2746 residuals+i,
2747 tree);
2748 }
2749
2750 // Initial targets
2751 n_targets = 3;
2752 for(i=0;i<n_targets;i++) targets[i] = tree->a_nodes[i]->b[0];
2753
2754 // Regraft each tip on the tree at most parsimonious position
2755 for(i=0;i<tree->n_otu-3;i++)
2756 {
2757 Set_Both_Sides(YES,tree);
2758 Lk(NULL,tree);
2759
2760 printf("\n. [%d/%d]",i,tree->n_otu-3);
2761
2762 tree->best_lnL = UNLIKELY;
2763 best_target = NULL;
2764 tg_idx = Permutate(n_targets);
2765
2766 for(j=0;j<n_targets;j++)
2767 {
2768 Graft_Subtree(targets[tg_idx[j]],
2769 tree->a_nodes[nd_idx[i]+3]->v[0],
2770 NULL,
2771 residuals[i],
2772 NULL,
2773 tree);
2774
2775 Update_PMat_At_Given_Edge(targets[tg_idx[j]],tree);
2776 Update_PMat_At_Given_Edge(tree->a_nodes[nd_idx[i]+3]->b[0],tree);
2777 Update_Partial_Lk(tree,residuals[i],tree->a_nodes[nd_idx[i]+3]->v[0]);
2778 Lk(residuals[i],tree);
2779
2780 if(tree->c_lnL > tree->best_lnL)
2781 {
2782 tree->best_lnL = tree->c_lnL;
2783 best_target = targets[tg_idx[j]];
2784 }
2785
2786 Prune_Subtree(tree->a_nodes[nd_idx[i]+3]->v[0],
2787 tree->a_nodes[nd_idx[i]+3],
2788 NULL,
2789 residuals+i,
2790 tree);
2791 }
2792
2793 assert(best_target);
2794
2795 Graft_Subtree(best_target,
2796 tree->a_nodes[nd_idx[i]+3]->v[0],
2797 NULL,
2798 residuals[i],
2799 NULL,
2800 tree);
2801
2802 n_opt = 0;
2803 do Optimize_Br_Len_Serie (2,tree); while(n_opt++ < 3);
2804
2805 targets[n_targets] = residuals[i];
2806 targets[n_targets+1] = tree->a_nodes[nd_idx[i]+3]->b[0];
2807
2808 Free(tg_idx);
2809 n_targets+=2;
2810 }
2811
2812 Round_Optimize(tree,5);
2813 PhyML_Fprintf(stderr,"\n. lk: %f",tree->c_lnL);
2814 Exit("\n");
2815
2816 Free(nd_idx);
2817 Free(residuals);
2818 Free(targets);
2819 }
2820
2821 //////////////////////////////////////////////////////////////
2822 //////////////////////////////////////////////////////////////
2823
2824
2825 /* |
2826 |
2827 |b
2828 |
2829 |d
2830 / \
2831 / \
2832 / \
2833 / \
2834 /v1 \v2
2835
2836 Set p_lk and sum_scale for subtrees with d, v1 and v2 as root,
2837 Pij for edges b, and the two edges connecting d to v1 and d to
2838 v2;
2839 Account for rooted trees.
2840 */
2841
Set_All_Partial_Lk(t_node ** n_v1,t_node ** n_v2,phydbl ** p_lk,int ** sum_scale,int ** p_lk_loc,phydbl ** Pij1,phydbl ** tPij1,phydbl ** p_lk_v1,int ** sum_scale_v1,phydbl ** Pij2,phydbl ** tPij2,phydbl ** p_lk_v2,int ** sum_scale_v2,t_node * d,t_edge * b,t_tree * tree,int * dest_p_idx,int * child1_p_idx,int * child2_p_idx,int * Pij1_idx,int * Pij2_idx)2842 void Set_All_Partial_Lk(t_node **n_v1, t_node **n_v2,
2843 phydbl **p_lk, int **sum_scale, int **p_lk_loc,
2844 phydbl **Pij1, phydbl **tPij1, phydbl **p_lk_v1, int **sum_scale_v1,
2845 phydbl **Pij2, phydbl **tPij2, phydbl **p_lk_v2, int **sum_scale_v2,
2846 t_node *d, t_edge *b, t_tree *tree
2847 #ifdef BEAGLE
2848 , int *dest_p_idx, int *child1_p_idx, int* child2_p_idx, int* Pij1_idx, int* Pij2_idx
2849 #endif
2850 )
2851 {
2852 unsigned int i;
2853
2854 assert(tree->is_mixt_tree == NO);
2855 assert(d->tax == NO);
2856
2857 if(tree->n_root == NULL || tree->ignore_root == YES)
2858 {
2859 /* Does d lie on the "left" or "right" of the branch? */
2860 if(d == b->left)
2861 {
2862 *p_lk = b->p_lk_left;
2863 *sum_scale = b->sum_scale_left;
2864 *p_lk_loc = b->p_lk_loc_left;
2865 #ifdef BEAGLE
2866 *dest_p_idx = b->p_lk_left_idx;
2867 #endif
2868 }
2869 else
2870 {
2871 *p_lk = b->p_lk_rght;
2872 *sum_scale = b->sum_scale_rght;
2873 *p_lk_loc = b->p_lk_loc_rght;
2874 #ifdef BEAGLE
2875 *dest_p_idx = b->p_lk_rght_idx;
2876 #endif
2877 }
2878
2879 *n_v1 = *n_v2 = NULL;
2880 for(i=0;i<3;++i)
2881 {
2882 if(d->b[i] != b)
2883 {
2884 if(!(*n_v1))
2885 {
2886 *n_v1 = d->v[i];
2887 #ifdef BEAGLE
2888 Set_Partial_Lk_One_Side(Pij1,tPij1,p_lk_v1,sum_scale_v1,d,d->b[i],tree,child1_p_idx,Pij1_idx);
2889 #else
2890 Set_Partial_Lk_One_Side(Pij1,tPij1,p_lk_v1,sum_scale_v1,d,d->b[i],tree);
2891 #endif
2892 }
2893 else if(!(*n_v2))
2894 {
2895 *n_v2 = d->v[i];
2896 #ifdef BEAGLE
2897 Set_Partial_Lk_One_Side(Pij2,tPij2,p_lk_v2,sum_scale_v2,d,d->b[i],tree,child2_p_idx,Pij2_idx);
2898 #else
2899 Set_Partial_Lk_One_Side(Pij2,tPij2,p_lk_v2,sum_scale_v2,d,d->b[i],tree);
2900 #endif
2901 }
2902 else assert(FALSE);
2903 }
2904 }
2905 }
2906 else
2907 {
2908 if(b == tree->e_root)
2909 {
2910 if(d == tree->n_root->v[1]) b = tree->n_root->b[1];
2911 else if(d == tree->n_root->v[2]) b = tree->n_root->b[2];
2912 else assert(FALSE);
2913 }
2914
2915 if(d == tree->n_root)
2916 {
2917 if(b == tree->n_root->b[1])
2918 {
2919 *p_lk = tree->n_root->b[1]->p_lk_left;
2920 *sum_scale = tree->n_root->b[1]->sum_scale_left;
2921 *p_lk_loc = tree->n_root->b[1]->p_lk_loc_left;
2922 #ifdef BEAGLE
2923 *dest_p_idx = tree->n_root->b[1]->p_lk_left_idx;
2924 #endif
2925 }
2926 else
2927 {
2928 *p_lk = tree->n_root->b[2]->p_lk_left;
2929 *sum_scale = tree->n_root->b[2]->sum_scale_left;
2930 *p_lk_loc = tree->n_root->b[2]->p_lk_loc_left;
2931 #ifdef BEAGLE
2932 *dest_p_idx = tree->n_root->b[2]->p_lk_left_idx;
2933 #endif
2934 }
2935
2936 *n_v1 = NULL;
2937 *Pij1 = NULL;
2938 *tPij1 = NULL;
2939 *p_lk_v1 = NULL;
2940 *sum_scale_v1 = NULL;
2941
2942 if(b == tree->n_root->b[1])
2943 {
2944 *n_v2 = tree->n_root->v[2];
2945 *Pij2 = tree->n_root->b[2]->Pij_rr;
2946 *tPij2 = tree->n_root->b[2]->tPij_rr;
2947 *p_lk_v2 = tree->n_root->b[2]->p_lk_rght;
2948 *sum_scale_v2 = tree->n_root->b[2]->sum_scale_rght;
2949 #ifdef BEAGLE
2950 *child2_p_idx = tree->n_root->b[2]->p_lk_rght_idx;
2951 *Pij2_idx = tree->n_root->b[2]->Pij_rr_idx;
2952 #endif
2953 }
2954 else if(b == tree->n_root->b[2])
2955 {
2956 *n_v2 = tree->n_root->v[1];
2957 *Pij2 = tree->n_root->b[1]->Pij_rr;
2958 *tPij2 = tree->n_root->b[1]->tPij_rr;
2959 *p_lk_v2 = tree->n_root->b[1]->p_lk_rght;
2960 *sum_scale_v2 = tree->n_root->b[1]->sum_scale_rght;
2961 #ifdef BEAGLE
2962 *child2_p_idx = tree->n_root->b[1]->p_lk_rght_idx;
2963 *Pij2_idx = tree->n_root->b[1]->Pij_rr_idx;
2964 #endif
2965 }
2966 else assert(FALSE);
2967 }
2968 else if(d == tree->n_root->v[1] || d == tree->n_root->v[2])
2969 {
2970 if(b == tree->n_root->b[1] || b == tree->n_root->b[2])
2971 {
2972 if(b == tree->n_root->b[1])
2973 {
2974 *p_lk = tree->n_root->b[1]->p_lk_rght;
2975 *sum_scale = tree->n_root->b[1]->sum_scale_rght;
2976 *p_lk_loc = tree->n_root->b[1]->p_lk_loc_left;
2977 #ifdef BEAGLE
2978 *dest_p_idx = tree->n_root->b[1]->p_lk_rght_idx;
2979 #endif
2980 }
2981 else
2982 {
2983 *p_lk = tree->n_root->b[2]->p_lk_rght;
2984 *sum_scale = tree->n_root->b[2]->sum_scale_rght;
2985 *p_lk_loc = tree->n_root->b[2]->p_lk_loc_rght;
2986 #ifdef BEAGLE
2987 *dest_p_idx = tree->n_root->b[2]->p_lk_rght_idx;
2988 #endif
2989 }
2990
2991 *n_v1 = *n_v2 = NULL;
2992 for(i=0;i<3;++i)
2993 {
2994 if(d->b[i] != tree->e_root)
2995 {
2996 if(!(*n_v1))
2997 {
2998 *n_v1 = d->v[i];
2999 #ifdef BEAGLE
3000 Set_Partial_Lk_One_Side(Pij1,tPij1,p_lk_v1,sum_scale_v1,d,d->b[i],tree,child1_p_idx,Pij1_idx);
3001 #else
3002 Set_Partial_Lk_One_Side(Pij1,tPij1,p_lk_v1,sum_scale_v1,d,d->b[i],tree);
3003 #endif
3004 }
3005 else
3006 {
3007 *n_v2 = d->v[i];
3008 #ifdef BEAGLE
3009 Set_Partial_Lk_One_Side(Pij2,tPij2,p_lk_v2,sum_scale_v2,d,d->b[i],tree,child2_p_idx,Pij2_idx);
3010 #else
3011 Set_Partial_Lk_One_Side(Pij2,tPij2,p_lk_v2,sum_scale_v2,d,d->b[i],tree);
3012 #endif
3013 }
3014 }
3015 }
3016 }
3017 else
3018 {
3019 if(d == b->left)
3020 {
3021 *p_lk = b->p_lk_left;
3022 *sum_scale = b->sum_scale_left;
3023 *p_lk_loc = b->p_lk_loc_left;
3024 #ifdef BEAGLE
3025 *dest_p_idx = b->p_lk_left_idx;
3026 #endif
3027 }
3028 else
3029 {
3030 *p_lk = b->p_lk_rght;
3031 *sum_scale = b->sum_scale_rght;
3032 *p_lk_loc = b->p_lk_loc_rght;
3033 #ifdef BEAGLE
3034 *dest_p_idx = b->p_lk_rght_idx;
3035 #endif
3036 }
3037
3038
3039 *n_v1 = tree->n_root;
3040 #ifdef BEAGLE
3041 Set_Partial_Lk_One_Side(Pij1,tPij1,p_lk_v1,sum_scale_v1,d,
3042 (d == tree->n_root->v[1])?
3043 (tree->n_root->b[1]):
3044 (tree->n_root->b[2]),
3045 tree,child1_p_idx,Pij1_idx);
3046 #else
3047 Set_Partial_Lk_One_Side(Pij1,tPij1,p_lk_v1,sum_scale_v1,d,
3048 (d == tree->n_root->v[1])?
3049 (tree->n_root->b[1]):
3050 (tree->n_root->b[2]),
3051 tree);
3052 #endif
3053 for(i=0;i<3;i++)
3054 {
3055 if(d->b[i] != tree->e_root && d->b[i] != b)
3056 {
3057 *n_v2 = d->v[i];
3058 #ifdef BEAGLE
3059 Set_Partial_Lk_One_Side(Pij2,tPij2,p_lk_v2,sum_scale_v2,d,d->b[i],tree,child2_p_idx,Pij2_idx);
3060 #else
3061 Set_Partial_Lk_One_Side(Pij2,tPij2,p_lk_v2,sum_scale_v2,d,d->b[i],tree);
3062 #endif
3063 break;
3064 }
3065 }
3066
3067 }
3068 }
3069 else
3070 {
3071 if(d == b->left)
3072 {
3073 *p_lk = b->p_lk_left;
3074 *sum_scale = b->sum_scale_left;
3075 *p_lk_loc = b->p_lk_loc_left;
3076 #ifdef BEAGLE
3077 *dest_p_idx = b->p_lk_left_idx;
3078 #endif
3079 }
3080 else
3081 {
3082 *p_lk = b->p_lk_rght;
3083 *sum_scale = b->sum_scale_rght;
3084 *p_lk_loc = b->p_lk_loc_rght;
3085 #ifdef BEAGLE
3086 *dest_p_idx = b->p_lk_rght_idx;
3087 #endif
3088 }
3089
3090 *n_v1 = *n_v2 = NULL;
3091 for(i=0;i<3;i++)
3092 {
3093 if(d->b[i] != b)
3094 {
3095 if(!(*n_v1))
3096 {
3097 *n_v1 = d->v[i];
3098 #ifdef BEAGLE
3099 Set_Partial_Lk_One_Side(Pij1,tPij1,p_lk_v1,sum_scale_v1,d,d->b[i],tree,child1_p_idx,Pij1_idx);
3100 #else
3101 Set_Partial_Lk_One_Side(Pij1,tPij1,p_lk_v1,sum_scale_v1,d,d->b[i],tree);
3102 #endif
3103 }
3104 else
3105 {
3106 *n_v2 = d->v[i];
3107 #ifdef BEAGLE
3108 Set_Partial_Lk_One_Side(Pij2,tPij2,p_lk_v2,sum_scale_v2,d,d->b[i],tree,child2_p_idx,Pij2_idx);
3109 #else
3110 Set_Partial_Lk_One_Side(Pij2,tPij2,p_lk_v2,sum_scale_v2,d,d->b[i],tree);
3111 #endif
3112 }
3113 }
3114 }
3115 }
3116 }
3117 }
3118
3119 //////////////////////////////////////////////////////////////
3120 //////////////////////////////////////////////////////////////
3121
3122 /* |
3123 |
3124 |d
3125 / \
3126 / \
3127 / \b
3128 / \
3129 / \x (either n_v1 or n_v2)
3130
3131 Returns p_lk and sum_scale for subtree with x as root, Pij for edge b
3132 */
3133
Set_Partial_Lk_One_Side(phydbl ** Pij,phydbl ** tPij,phydbl ** p_lk,int ** sum_scale,t_node * d,t_edge * b,t_tree * tree,int * child_p_idx,int * Pij_idx)3134 void Set_Partial_Lk_One_Side(phydbl **Pij, phydbl **tPij, phydbl **p_lk, int **sum_scale, t_node *d, t_edge *b, t_tree *tree
3135 #ifdef BEAGLE
3136 , int* child_p_idx, int* Pij_idx
3137 #endif
3138 )
3139 {
3140
3141 if(Pij != NULL)
3142 {
3143 *Pij = b->Pij_rr;
3144 *tPij = b->tPij_rr;
3145 #ifdef BEAGLE
3146 *Pij_idx = b->Pij_rr_idx;
3147 #endif
3148 }
3149
3150 if(d->tax == NO)
3151 {
3152 if(d == b->left) // if d is on the left of b, then d's neighbor is on the right
3153 {
3154 *p_lk = (b->rght->tax == YES) ? b->p_lk_tip_r : b->p_lk_rght;
3155 *sum_scale = b->sum_scale_rght;
3156 #ifdef BEAGLE
3157 *child_p_idx = b->rght->tax? b->p_lk_tip_idx: b->p_lk_rght_idx;
3158 #endif
3159 assert(*p_lk);
3160 }
3161 else
3162 {
3163 *p_lk = b->p_lk_left;
3164 *sum_scale = b->sum_scale_left;
3165 #ifdef BEAGLE
3166 *child_p_idx = b->rght->tax? b->p_lk_tip_idx: b->p_lk_left_idx;
3167 #endif
3168 assert(*p_lk);
3169 }
3170 }
3171 else
3172 {
3173 #ifdef BEAGLE
3174 Warn_And_Exit(TODO_BEAGLE);
3175 #endif
3176 *p_lk = NULL;
3177 *sum_scale = NULL;
3178 }
3179 }
3180
3181 //////////////////////////////////////////////////////////////
3182 //////////////////////////////////////////////////////////////
3183
Switch_Partial_Lk_Post(t_node * a,t_node * d,t_edge * b,short int yesno,t_tree * tree)3184 void Switch_Partial_Lk_Post(t_node *a, t_node *d, t_edge *b, short int yesno, t_tree *tree)
3185 {
3186 if(a == tree->n_root) assert(FALSE);
3187 if(d->tax == NO)
3188 {
3189 int i;
3190
3191 for(i=0;i<3;++i)
3192 if(d->v[i] != a)
3193 Switch_Partial_Lk_Post(d,d->v[i],d->b[i],yesno,tree);
3194 }
3195
3196 if(b->left == d) b->update_partial_lk_left = yesno;
3197 else if(b->rght == d) b->update_partial_lk_rght = yesno;
3198 else assert(FALSE);
3199
3200 return;
3201 }
3202
3203
3204 //////////////////////////////////////////////////////////////
3205 //////////////////////////////////////////////////////////////
3206
Switch_Partial_Lk_Pre(t_node * a,t_node * d,t_edge * b,short int yesno,t_tree * tree)3207 void Switch_Partial_Lk_Pre(t_node *a, t_node *d, t_edge *b, short int yesno, t_tree *tree)
3208 {
3209 if(a == tree->n_root) assert(FALSE);
3210 if(d->tax == YES) return;
3211 else
3212 {
3213 int i;
3214
3215 for(i=0;i<3;++i)
3216 {
3217 if(d->v[i] != a)
3218 {
3219 if(d->b[i]->left == d) d->b[i]->update_partial_lk_left = yesno;
3220 else if(d->b[i]->rght == d) d->b[i]->update_partial_lk_rght = yesno;
3221 else assert(FALSE);
3222 }
3223 }
3224
3225 for(i=0;i<3;++i)
3226 if(d->v[i] != a)
3227 Switch_Partial_Lk_Pre(d,d->v[i],d->b[i],yesno,tree);
3228 }
3229
3230 return;
3231 }
3232
3233 //////////////////////////////////////////////////////////////
3234 //////////////////////////////////////////////////////////////
3235
Partial_Lk_Inin(const phydbl * Pij1,const phydbl * plk1,const phydbl * Pij2,const phydbl * plk2,const int ns,phydbl * plk0)3236 void Partial_Lk_Inin(const phydbl *Pij1, const phydbl *plk1, const phydbl *Pij2, const phydbl *plk2, const int ns, phydbl *plk0)
3237 {
3238 unsigned int i,j;
3239
3240
3241 for(i=0;i<ns;++i) if(plk1[i] > 1.0 || plk1[i] < 1.0 || plk2[i] > 1.0 || plk2[i] < 1.0) break;
3242
3243 if(i != ns)
3244 {
3245 for(i=0;i<ns;++i)
3246 {
3247 phydbl u1 = 0.0;
3248 phydbl u2 = 0.0;
3249
3250 for(j=0;j<ns;++j)
3251 {
3252 u1 += Pij1[j] * plk1[j];
3253 u2 += Pij2[j] * plk2[j];
3254 }
3255
3256 Pij1 += ns;
3257 Pij2 += ns;
3258 plk0[i] = u1*u2;
3259 }
3260 }
3261 else
3262 {
3263 for(i=0;i<ns;++i) plk0[i] = 1.0;
3264 }
3265 }
3266
3267 //////////////////////////////////////////////////////////////
3268 //////////////////////////////////////////////////////////////
3269
Partial_Lk_Exex(const phydbl * Pij1,const int state1,const phydbl * Pij2,const int state2,const int ns,phydbl * plk0)3270 void Partial_Lk_Exex(const phydbl *Pij1, const int state1, const phydbl *Pij2, const int state2, const int ns, phydbl *plk0)
3271 {
3272 unsigned int i;
3273 for(i=0;i<ns;++i)
3274 {
3275 plk0[i] = Pij1[state1]*Pij2[state2];
3276 Pij1 += ns;
3277 Pij2 += ns;
3278 }
3279 }
3280
3281 //////////////////////////////////////////////////////////////
3282 //////////////////////////////////////////////////////////////
3283
Partial_Lk_Exin(const phydbl * Pij1,const int state1,const phydbl * Pij2,const phydbl * plk2,const int ns,phydbl * plk0)3284 void Partial_Lk_Exin(const phydbl *Pij1, const int state1, const phydbl *Pij2, const phydbl *plk2, const int ns, phydbl *plk0)
3285 {
3286 unsigned int i,j;
3287
3288 for(i=0;i<ns;++i)
3289 {
3290 phydbl u2 = 0.0;
3291 for(j=0;j<ns;++j) u2 += Pij2[j] * plk2[j];
3292 plk0[i] = Pij1[state1]*u2;
3293 Pij1 += ns;
3294 Pij2 += ns;
3295 }
3296 }
3297
3298 //////////////////////////////////////////////////////////////
3299 //////////////////////////////////////////////////////////////
3300