1 /*
2
3 PHYML : a program that computes maximum likelihood phylogenies from
4 DNA or AA homologous sequences
5
6 Copyright (C) Stephane Guindon. Oct 2003 onward
7
8 All parts of the source except where indicated are distributed under
9 the GNU public licence. See http://www.opensource.org for details.
10
11 */
12
13 /*
14
15 Implementation of aLRT branch tests, and 5-branchs NNI research optimization.
16
17 Authors : Jean-Francois Dufayard & Stephane Guindon.
18
19 */
20
21 #include "alrt.h"
22
23 //////////////////////////////////////////////////////////////
24 //////////////////////////////////////////////////////////////
25
26 /*
27 * Check every testable branch of the tree,
28 * check for NNIs, optimizing 5 branches,
29 * param tree : the tree to check
30 */
31
Check_NNI_Five_Branches(t_tree * tree)32 int Check_NNI_Five_Branches(t_tree *tree)
33 {
34 int best_edge;
35 phydbl best_gain;
36 int best_config;
37 int i;
38 int better_found; /* = 1 if a phylogeny with greater likelihood than current one was found */
39 int result;
40 phydbl init_lnL;
41
42 init_lnL = UNLIKELY;
43 better_found = YES;
44
45 //While there is at least one NNI to do...
46 while(better_found == YES)
47 {
48 Update_Dirs(tree);
49
50 //Interface output
51 if(tree->verbose > VL0 && tree->io->quiet == NO) PhyML_Printf("\n\n. Checking for NNIs, optimizing five branches...\n");
52
53 better_found = 0;
54 result = -1;
55 best_edge = -1;
56 best_gain = tree->mod->s_opt->min_diff_lk_local;
57 best_config = -1;
58
59 MIXT_Set_Alias_Subpatt(YES,tree);
60 Set_Both_Sides(YES,tree);
61 init_lnL = Lk(NULL,tree);
62 MIXT_Set_Alias_Subpatt(NO,tree);
63
64
65 //For every branch
66 for(i=0;i<2*tree->n_otu-3;++i)
67 {
68 //if this branch is not terminal
69 if((!tree->a_edges[i]->left->tax) && (!tree->a_edges[i]->rght->tax))
70 {
71 result = -1;
72
73 //Try every NNIs for the tested branch
74 result = NNI_Neigh_BL(tree->a_edges[i],tree);
75
76 //Look for possible NNI to do, and check if it is the best one
77 switch(result)
78 {
79 case 1 : /* lk1 > lk0 > lk2 */
80 {
81 if((tree->a_edges[i]->nni->lk0 - tree->a_edges[i]->nni->lk1) < best_gain)
82 {
83 better_found = YES;
84 best_edge = i;
85 best_gain = tree->a_edges[i]->nni->lk0-tree->a_edges[i]->nni->lk1;
86 best_config = 1;
87 }
88 break;
89 }
90 case 2 : /* lk2 > lk0 > lk1 */
91 {
92 if((tree->a_edges[i]->nni->lk0 - tree->a_edges[i]->nni->lk2) < best_gain)
93 {
94 better_found = YES;
95 best_edge = i;
96 best_gain = tree->a_edges[i]->nni->lk0-tree->a_edges[i]->nni->lk2;
97 best_config = 2;
98 }
99 break;
100 }
101 case 3 : /* lk1 > lk2 > lk0 */
102 {
103 if((tree->a_edges[i]->nni->lk2 - tree->a_edges[i]->nni->lk1) < best_gain)
104 {
105 better_found = YES;
106 best_edge = i;
107 best_gain = tree->a_edges[i]->nni->lk2-tree->a_edges[i]->nni->lk1;
108 best_config = 1;
109 }
110 break;
111 }
112 case 4 : /* lk2 > lk1 > lk0 */
113 {
114 if((tree->a_edges[i]->nni->lk1 - tree->a_edges[i]->nni->lk2) < best_gain)
115 {
116 better_found = YES;
117 best_edge = i;
118 best_gain = tree->a_edges[i]->nni->lk1-tree->a_edges[i]->nni->lk2;
119 best_config = 2;
120 }
121 break;
122 }
123 default : /* lk2 = lk1 = lk0 */
124 {
125 if(best_gain > .0) best_gain = .0;
126 break;
127 }
128 }
129 }
130 }
131
132
133 if((tree->c_lnL < init_lnL - tree->mod->s_opt->min_diff_lk_local) || (tree->c_lnL > init_lnL + tree->mod->s_opt->min_diff_lk_local))
134 {
135 PhyML_Fprintf(stderr,"\n\n== tree->c_lnL = %f init_lnL = %f.",tree->c_lnL,init_lnL);
136 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d\n\n.\n",__FILE__,__LINE__);
137 Warn_And_Exit("\n");
138 }
139
140 //Don't do any NNI if the user doesn't want to optimize topology
141 if(!tree->mod->s_opt->opt_topo) better_found = NO;
142 /* if(FABS(best_gain) <= tree->mod->s_opt->min_diff_lk_move) better_found = 0; */
143
144 //If there is one swap to do, make the best one.
145 if(better_found)
146 {
147 Make_Target_Swap(tree,tree->a_edges[best_edge],best_config);
148
149 MIXT_Set_Alias_Subpatt(YES,tree);
150 Lk(NULL,tree);
151 MIXT_Set_Alias_Subpatt(YES,tree);
152
153 if(tree->c_lnL < init_lnL)
154 {
155 PhyML_Fprintf(stderr,"\n\n== tree->c_lnL = %f init_lnL = %f.",tree->c_lnL,init_lnL);
156 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d\n\n.\n",__FILE__,__LINE__);
157 Exit("\n");
158 }
159
160 if(tree->verbose > VL0 && tree->io->quiet == NO) Print_Lk(tree,"[Topology ]");
161
162 if(FABS(tree->c_lnL - init_lnL) < tree->mod->s_opt->min_diff_lk_move) return 0;
163 return 1;
164 }
165 }
166 return 0;
167 }
168
169 //////////////////////////////////////////////////////////////
170 //////////////////////////////////////////////////////////////
171
172 /* Compute aLRT supports */
aLRT(t_tree * tree)173 void aLRT(t_tree *tree)
174 {
175 int i;
176 char *method;
177
178 // Print_All_Edge_Likelihoods(tree);
179 Unscale_Br_Len_Multiplier_Tree(tree);
180 Br_Len_Not_Involving_Invar(tree);
181 // Print_All_Edge_Likelihoods(tree);
182 /* aLRT support will label each internal branch */
183 tree->io->print_support_val = YES;
184
185 /* The topology will not be modified when assessing the branch support. We make sure that it will
186 not be modified afterwards by locking the topology */
187
188 method = (char *)mCalloc(100,sizeof(char));
189
190 switch(tree->io->ratio_test)
191 {
192 case ALRTCHI2: { strcpy(method,"aLRT"); break; }
193 case MINALRTCHI2SH: { strcpy(method,"aLRT"); break; }
194 case ALRTSTAT: { strcpy(method,"aLRT"); break; }
195 case SH: { strcpy(method,"SH"); break; }
196 case ABAYES: { strcpy(method,"aBayes"); break; }
197 default : return;
198 }
199
200 if(tree->io->quiet == NO) PhyML_Printf("\n\n. Calculating fast branch supports (using '%s').",method);
201 Free(method);
202
203 MIXT_Set_Alias_Subpatt(YES,tree);
204 Set_Both_Sides(YES,tree);
205 // Print_All_Edge_Likelihoods(tree);
206 Lk(NULL,tree);
207 // Print_All_Edge_Likelihoods(tree);
208 MIXT_Set_Alias_Subpatt(NO,tree);
209 Update_Dirs(tree);
210
211 for(i=0;i<2*tree->n_otu-3;++i)
212 {
213 if((!tree->a_edges[i]->left->tax) && (!tree->a_edges[i]->rght->tax))
214 {
215 /* Compute likelihoods for each of the three configuration */
216 NNI_Neigh_BL(tree->a_edges[i],tree);
217 /* Compute the corresponding statistical support */
218 Compute_Likelihood_Ratio_Test(tree->a_edges[i],tree);
219 }
220 }
221
222 tree->lock_topo = YES;
223
224 Br_Len_Involving_Invar(tree);
225 Rescale_Br_Len_Multiplier_Tree(tree);
226
227 }
228
229 //////////////////////////////////////////////////////////////
230 //////////////////////////////////////////////////////////////
231
232 /*
233 * Launch one branch testing,
234 * analyse the result
235 * and convert supports as wished by the user.
236 * param tree : the tree to check
237 * param tested_t_edge : the tested t_edge of the tree
238 * param old_loglk : the initial likelihood, before any aLRT analysis
239 * param isBoot : 1 if used from the Bootstrap procedure, 0 if not
240 * return an integer, informative to analyse the results and potential NNIs to do
241 */
Compute_Likelihood_Ratio_Test(t_edge * tested_edge,t_tree * tree)242 int Compute_Likelihood_Ratio_Test(t_edge *tested_edge, t_tree *tree)
243 {
244 int result=0;
245
246 tested_edge->ratio_test = 0.0;
247 tested_edge->alrt_statistic = -1.0;
248
249 if((tested_edge->nni->lk0 > tested_edge->nni->lk1) && (tested_edge->nni->lk0 > tested_edge->nni->lk2))
250 {
251 if(tested_edge->nni->lk1 < tested_edge->nni->lk2)
252 {
253 //lk0 > lk2 > lk1
254 tested_edge->alrt_statistic = 2*(tested_edge->nni->lk0 - tested_edge->nni->lk2);
255 }
256 else
257 {
258 //lk0 > lk1 >= lk2
259 tested_edge->alrt_statistic = 2*(tested_edge->nni->lk0 - tested_edge->nni->lk1);
260 }
261
262 if (tested_edge->alrt_statistic < 0.0)
263 {
264 tested_edge->alrt_statistic = 0.0;
265 tested_edge->ratio_test = 0.0;
266 }
267 else
268 {
269 //aLRT statistic is valid, compute the branch support
270 if (tree->io->ratio_test == ALRTCHI2)
271 {
272 tested_edge->ratio_test = Statistics_To_Probabilities(tested_edge->alrt_statistic);
273 }
274 else if(tree->io->ratio_test == MINALRTCHI2SH)
275 {
276 phydbl sh_support;
277 phydbl param_support;
278
279 sh_support = Statistics_To_SH(tree);
280 param_support = Statistics_To_Probabilities(tested_edge->alrt_statistic);
281
282 if(sh_support < param_support) tested_edge->ratio_test = sh_support;
283 else tested_edge->ratio_test = param_support;
284 }
285
286 else if(tree->io->ratio_test == ALRTSTAT)
287 {
288 tested_edge->ratio_test=tested_edge->alrt_statistic;
289 }
290 else if(tree->io->ratio_test == SH)
291 {
292 tested_edge->ratio_test = Statistics_To_SH(tree);
293 }
294 else if(tree->io->ratio_test == ABAYES)
295 {
296 phydbl Kp0,Kp1,Kp2,logK;
297
298 logK = 1. - MAX(MAX(tested_edge->nni->lk0,tested_edge->nni->lk1),tested_edge->nni->lk2);
299
300 Kp0 = exp(tested_edge->nni->lk0+logK);
301 Kp1 = exp(tested_edge->nni->lk1+logK);
302 Kp2 = exp(tested_edge->nni->lk2+logK);
303 tested_edge->ratio_test = Kp0/(Kp0+Kp1+Kp2);
304 }
305 }
306 }
307 //statistic is not valid, give the negative statistics if wished, or a zero support.
308 else if((tested_edge->nni->lk1 > tested_edge->nni->lk0) && (tested_edge->nni->lk1 > tested_edge->nni->lk2))
309 {
310 /* tested_edge->ratio_test = 2*(tested_edge->nni->lk0-tested_edge->nni->lk1); */
311 tested_edge->ratio_test = 0.0;
312 if(tree->io->ratio_test > 1) tested_edge->alrt_statistic = 0.0;
313 }
314 else if((tested_edge->nni->lk2 > tested_edge->nni->lk0) && (tested_edge->nni->lk2 > tested_edge->nni->lk1))
315 {
316 /* tested_edge->ratio_test = 2*(tested_edge->nni->lk0-tested_edge->nni->lk2); */
317 tested_edge->ratio_test = 0.0;
318 if(tree->io->ratio_test > 1) tested_edge->alrt_statistic = 0.0;
319 }
320 else // lk0 ~ lk1 ~ lk2
321 {
322 tested_edge->ratio_test = 0.0;
323 if(tree->io->ratio_test > 1) tested_edge->ratio_test = 0.0;
324 }
325 return result;
326 }
327
328 //////////////////////////////////////////////////////////////
329 //////////////////////////////////////////////////////////////
330
331 /*
332 * Test the 3 NNI positions for one branch.
333 * param tree : the tree to check
334 * param tested_t_edge : the tested t_edge of the tree
335 * param old_loglk : the initial likelihood, before any aLRT analysis
336 * param isBoot : 1 if used from the Bootstrap procedure, 0 if not
337 */
NNI_Neigh_BL(t_edge * b_fcus,t_tree * tree)338 int NNI_Neigh_BL(t_edge *b_fcus, t_tree *tree)
339 {
340 /*!
341 syntax : (node) [edge]
342 (left_1) . .(right_1)
343 \ (left) (right) /
344 \._____________./
345 / [b_fcus] \
346 / \
347 (left_2) . .(right_2)
348
349 */
350
351 int site;
352 t_node *v1,*v2,*v3,*v4;
353 t_edge *e1,*e2,*e3,*e4;
354 scalar_dbl *len_e1,*len_e2,*len_e3,*len_e4;
355 scalar_dbl *var_e1,*var_e2,*var_e3,*var_e4;
356 phydbl lk0, lk1, lk2;
357 scalar_dbl *l_init,*v_init;
358 phydbl lk_init, lk_temp;
359 int i;
360 int result;
361 void *buff;
362
363 result = 0;
364
365
366 /*! Initialization */
367 l_init = Duplicate_Scalar_Dbl(b_fcus->l);
368 v_init = Duplicate_Scalar_Dbl(b_fcus->l_var);
369 lk_init = tree->c_lnL;
370 lk_temp = UNLIKELY;
371 lk0 = lk1 = lk2 = UNLIKELY;
372 v1 = v2 = v3 = v4 = NULL;
373 Init_NNI_Score(0.0,b_fcus,tree);
374
375 /*! vertices */
376 v1 = b_fcus->left->v[b_fcus->l_v1];
377 v2 = b_fcus->left->v[b_fcus->l_v2];
378 v3 = b_fcus->rght->v[b_fcus->r_v1];
379 v4 = b_fcus->rght->v[b_fcus->r_v2];
380
381 if(v1->num < v2->num)
382 {
383 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
384 Exit("");
385 }
386
387 if(v3->num < v4->num)
388 {
389 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
390 Exit("");
391 }
392
393 /*! edges */
394 e1 = b_fcus->left->b[b_fcus->l_v1];
395 e2 = b_fcus->left->b[b_fcus->l_v2];
396 e3 = b_fcus->rght->b[b_fcus->r_v1];
397 e4 = b_fcus->rght->b[b_fcus->r_v2];
398
399 /*! record initial branch lengths */
400 len_e1 = Duplicate_Scalar_Dbl(e1->l);
401 len_e2 = Duplicate_Scalar_Dbl(e2->l);
402 len_e3 = Duplicate_Scalar_Dbl(e3->l);
403 len_e4 = Duplicate_Scalar_Dbl(e4->l);
404 var_e1 = Duplicate_Scalar_Dbl(e1->l_var);
405 var_e2 = Duplicate_Scalar_Dbl(e2->l_var);
406 var_e3 = Duplicate_Scalar_Dbl(e3->l_var);
407 var_e4 = Duplicate_Scalar_Dbl(e4->l_var);
408
409 /*! Optimize branch lengths and update likelihoods for
410 the original configuration.
411 */
412 do
413 {
414 lk0 = lk_temp;
415
416 for(i=0;i<3;i++) //a branch has three nodes
417 if(b_fcus->left->v[i] != b_fcus->rght) //Only consider left_1 and left_2
418 {
419 Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
420 lk_temp = Br_Len_Opt(&(b_fcus->left->b[i]->l->v),b_fcus->left->b[i],tree);
421 }
422
423 Update_Partial_Lk(tree,b_fcus,b_fcus->left);
424 lk_temp = Br_Len_Opt(&(b_fcus->l->v),b_fcus,tree);
425
426 for(i=0;i<3;i++)
427 if(b_fcus->rght->v[i] != b_fcus->left) //Only consider right_1 and right_2
428 {
429 Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
430 lk_temp = Br_Len_Opt(&(b_fcus->rght->b[i]->l->v),b_fcus->rght->b[i],tree);
431 }
432 Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
433
434 if(lk_temp < lk0 - tree->mod->s_opt->min_diff_lk_local)
435 {
436 PhyML_Fprintf(stderr,"\n. lk_temp = %f lk0 = %f\n",lk_temp,lk0);
437 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
438 Exit("\n");
439 }
440 }
441 while(fabs(lk_temp-lk0) > tree->mod->s_opt->min_diff_lk_global);
442 //until no significative improvement is detected
443
444 lk0 = tree->c_lnL;
445
446 buff = (t_tree *)tree;
447 do
448 {
449 for(site=0;site<tree->n_pattern;site++) tree->log_lks_aLRT[0][site] = tree->c_lnL_sorted[site];
450 tree = tree->next_mixt;
451 }
452 while(tree);
453
454 tree = (t_tree *)buff;
455
456 /*! Go back to initial branch lengths */
457 Copy_Scalar_Dbl(len_e1,e1->l);
458 Copy_Scalar_Dbl(len_e2,e2->l);
459 Copy_Scalar_Dbl(len_e3,e3->l);
460 Copy_Scalar_Dbl(len_e4,e4->l);
461 Copy_Scalar_Dbl(l_init,b_fcus->l);
462 Copy_Scalar_Dbl(var_e1,e1->l_var);
463 Copy_Scalar_Dbl(var_e2,e2->l_var);
464 Copy_Scalar_Dbl(var_e3,e3->l_var);
465 Copy_Scalar_Dbl(var_e4,e4->l_var);
466 Copy_Scalar_Dbl(v_init,b_fcus->l_var);
467
468 Update_PMat_At_Given_Edge(e1,tree);
469 Update_PMat_At_Given_Edge(e2,tree);
470 Update_PMat_At_Given_Edge(e3,tree);
471 Update_PMat_At_Given_Edge(e4,tree);
472 Update_PMat_At_Given_Edge(b_fcus,tree);
473
474 /* Sanity check */
475 MIXT_Set_Alias_Subpatt(YES,tree);
476 Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
477 Update_Partial_Lk(tree,b_fcus,b_fcus->left);
478 lk_temp = Lk(b_fcus,tree);
479 MIXT_Set_Alias_Subpatt(NO,tree);
480 if((lk_temp > lk_init + tree->mod->s_opt->min_diff_lk_local) || (lk_temp < lk_init - tree->mod->s_opt->min_diff_lk_local))
481 {
482 PhyML_Fprintf(stderr,"\n. lk_temp = %f lk_init = %f",lk_temp,lk_init);
483 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
484 Exit("\n");
485 }
486
487
488 /*! Do first possible swap */
489 Swap(v2,b_fcus->left,b_fcus->rght,v3,tree);
490 MIXT_Set_Alias_Subpatt(YES,tree);
491 Set_Both_Sides(YES,tree);
492 Update_Partial_Lk(tree,b_fcus,b_fcus->left);
493 Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
494 lk1 = Lk(b_fcus,tree);
495 MIXT_Set_Alias_Subpatt(NO,tree);
496
497 for(i=0;i<3;i++)
498 if(b_fcus->left->v[i] != b_fcus->rght)
499 Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
500
501 for(i=0;i<3;i++)
502 if(b_fcus->rght->v[i] != b_fcus->left)
503 Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
504
505
506 /*! Optimize branch lengths and update likelihoods */
507 lk_temp = UNLIKELY;
508 do
509 {
510 lk1 = lk_temp;
511
512 for(i=0;i<3;i++)
513 if(b_fcus->left->v[i] != b_fcus->rght)
514 {
515 Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
516 lk_temp = Br_Len_Opt(&(b_fcus->left->b[i]->l->v),b_fcus->left->b[i],tree);
517 }
518
519 Update_Partial_Lk(tree,b_fcus,b_fcus->left);
520 lk_temp = Br_Len_Opt(&(b_fcus->l->v),b_fcus,tree);
521
522 for(i=0;i<3;i++)
523 if(b_fcus->rght->v[i] != b_fcus->left)
524 {
525 Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
526 lk_temp = Br_Len_Opt(&(b_fcus->rght->b[i]->l->v),b_fcus->rght->b[i],tree);
527 }
528
529 Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
530
531 if(lk_temp < lk1 - tree->mod->s_opt->min_diff_lk_local)
532 {
533 PhyML_Fprintf(stderr,"\n. lk_temp = %f lk1 = %f\n",lk_temp,lk1);
534 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
535 Exit("");
536 }
537 }
538 while(FABS(lk_temp-lk1) > tree->mod->s_opt->min_diff_lk_global);
539 /*! until no significative improvement is detected */
540
541 lk1 = tree->c_lnL;
542
543 /*! record likelihood of each site, in order to compute branch supports */
544
545 buff = (t_tree *)tree;
546 do
547 {
548 for(site=0;site<tree->n_pattern;site++) tree->log_lks_aLRT[1][site]= tree->c_lnL_sorted[site];
549 tree = tree->next_mixt;
550 }
551 while(tree);
552 tree = (t_tree *)buff;
553
554 Swap(v3,b_fcus->left,b_fcus->rght,v2,tree);
555
556 /*! Go back to initial branch lengths */
557 Copy_Scalar_Dbl(len_e1,e1->l);
558 Copy_Scalar_Dbl(len_e2,e2->l);
559 Copy_Scalar_Dbl(len_e3,e3->l);
560 Copy_Scalar_Dbl(len_e4,e4->l);
561 Copy_Scalar_Dbl(l_init,b_fcus->l);
562 Copy_Scalar_Dbl(var_e1,e1->l_var);
563 Copy_Scalar_Dbl(var_e2,e2->l_var);
564 Copy_Scalar_Dbl(var_e3,e3->l_var);
565 Copy_Scalar_Dbl(var_e4,e4->l_var);
566 Copy_Scalar_Dbl(v_init,b_fcus->l_var);
567
568 Update_PMat_At_Given_Edge(e1,tree);
569 Update_PMat_At_Given_Edge(e2,tree);
570 Update_PMat_At_Given_Edge(e3,tree);
571 Update_PMat_At_Given_Edge(e4,tree);
572 Update_PMat_At_Given_Edge(b_fcus,tree);
573
574
575 /* Sanity check */
576 MIXT_Set_Alias_Subpatt(YES,tree);
577 Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
578 Update_Partial_Lk(tree,b_fcus,b_fcus->left);
579 lk_temp = Lk(b_fcus,tree);
580 MIXT_Set_Alias_Subpatt(NO,tree);
581 if((lk_temp > lk_init + tree->mod->s_opt->min_diff_lk_local) || (lk_temp < lk_init - tree->mod->s_opt->min_diff_lk_local))
582 {
583 PhyML_Fprintf(stderr,"\n== lk_temp = %f lk_init = %f",lk_temp,lk_init);
584 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
585 Warn_And_Exit("");
586 }
587
588 /*! do the second possible swap */
589 Swap(v2,b_fcus->left,b_fcus->rght,v4,tree);
590 Set_Both_Sides(YES,tree);
591 MIXT_Set_Alias_Subpatt(YES,tree);
592 Update_Partial_Lk(tree,b_fcus,b_fcus->left);
593 Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
594 lk2 = Lk(b_fcus,tree);
595 MIXT_Set_Alias_Subpatt(NO,tree);
596
597 for(i=0;i<3;i++)
598 if(b_fcus->left->v[i] != b_fcus->rght)
599 Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
600
601 for(i=0;i<3;i++)
602 if(b_fcus->rght->v[i] != b_fcus->left)
603 Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
604
605 /*! Optimize branch lengths and update likelihoods */
606 lk_temp = UNLIKELY;
607 do
608 {
609 lk2 = lk_temp;
610
611 for(i=0;i<3;i++)
612 if(b_fcus->left->v[i] != b_fcus->rght)
613 {
614 Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
615 lk_temp = Br_Len_Opt(&(b_fcus->left->b[i]->l->v),b_fcus->left->b[i],tree);
616
617 if(lk_temp < lk2 - tree->mod->s_opt->min_diff_lk_local)
618 {
619 PhyML_Fprintf(stderr,"\n== l: %f var:%f",b_fcus->left->b[i]->l->v,b_fcus->left->b[i]->l_var->v);
620 PhyML_Fprintf(stderr,"\n== lk_temp = %f lk2 = %f",lk_temp,lk2);
621 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d",__FILE__,__LINE__);
622 Exit("\n");
623 }
624 }
625 Update_Partial_Lk(tree,b_fcus,b_fcus->left);
626 lk_temp = Br_Len_Opt(&(b_fcus->l->v),b_fcus,tree);
627
628 if(lk_temp < lk2 - tree->mod->s_opt->min_diff_lk_local)
629 {
630 PhyML_Fprintf(stderr,"\n== l: %f var:%f",b_fcus->l->v,b_fcus->l_var->v);
631 PhyML_Fprintf(stderr,"\n== lk_temp = %f lk2 = %f",lk_temp,lk2);
632 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d",__FILE__,__LINE__);
633 Exit("\n");
634 }
635
636 for(i=0;i<3;i++)
637 if(b_fcus->rght->v[i] != b_fcus->left)
638 {
639 Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
640 lk_temp = Br_Len_Opt(&(b_fcus->rght->b[i]->l->v),b_fcus->rght->b[i],tree);
641
642 if(lk_temp < lk2 - tree->mod->s_opt->min_diff_lk_local)
643 {
644 Set_Both_Sides(YES,tree);
645 Lk(b_fcus,tree);
646 Check_Lk_At_Given_Edge(YES,tree);
647 PhyML_Fprintf(stderr,"\n== l: %f var:%f",b_fcus->rght->b[i]->l->v,b_fcus->rght->b[i]->l_var->v);
648 PhyML_Fprintf(stderr,"\n== lk_temp = %f lk2 = %f",lk_temp,lk2);
649 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d",__FILE__,__LINE__);
650 Exit("\n");
651 }
652 }
653
654 Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
655
656 }
657 while(FABS(lk_temp-lk2) > tree->mod->s_opt->min_diff_lk_global);
658 //until no significative improvement is detected
659
660 lk2 = tree->c_lnL;
661
662 /*! save likelihood of each sites, in order to compute branch supports */
663 buff = (t_tree *)tree;
664 do
665 {
666 for(site=0;site<tree->n_pattern;site++) tree->log_lks_aLRT[2][site]= tree->c_lnL_sorted[site];
667 tree = tree->next_mixt;
668 }
669 while(tree);
670 tree = (t_tree *)buff;
671
672
673 /*! undo the swap */
674 Swap(v4,b_fcus->left,b_fcus->rght,v2,tree);
675 /***********/
676
677 /*! restore the initial branch length values */
678 Copy_Scalar_Dbl(len_e1,e1->l);
679 Copy_Scalar_Dbl(len_e2,e2->l);
680 Copy_Scalar_Dbl(len_e3,e3->l);
681 Copy_Scalar_Dbl(len_e4,e4->l);
682 Copy_Scalar_Dbl(l_init,b_fcus->l);
683 Copy_Scalar_Dbl(var_e1,e1->l_var);
684 Copy_Scalar_Dbl(var_e2,e2->l_var);
685 Copy_Scalar_Dbl(var_e3,e3->l_var);
686 Copy_Scalar_Dbl(var_e4,e4->l_var);
687 Copy_Scalar_Dbl(v_init,b_fcus->l_var);
688
689
690 /*! recompute likelihoods */
691 Update_PMat_At_Given_Edge(e1,tree);
692 Update_PMat_At_Given_Edge(e2,tree);
693 Update_PMat_At_Given_Edge(e3,tree);
694 Update_PMat_At_Given_Edge(e4,tree);
695 Update_PMat_At_Given_Edge(b_fcus,tree);
696
697 MIXT_Set_Alias_Subpatt(YES,tree);
698 Update_Partial_Lk(tree,b_fcus,b_fcus->left);
699 Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
700
701 for(i=0;i<3;i++)
702 if(b_fcus->left->v[i] != b_fcus->rght)
703 Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
704
705 for(i=0;i<3;i++)
706 if(b_fcus->rght->v[i] != b_fcus->left)
707 Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
708
709 MIXT_Set_Alias_Subpatt(NO,tree);
710
711 lk_temp = Lk(b_fcus,tree);
712
713 if((lk_temp > lk_init + tree->mod->s_opt->min_diff_lk_local) || (lk_temp < lk_init - tree->mod->s_opt->min_diff_lk_local))
714 {
715 PhyML_Fprintf(stderr,"\n== lk_temp = %f lk_init = %f",lk_temp,lk_init);
716 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
717 Warn_And_Exit("");
718 }
719
720
721 //save likelihoods in NNI structures
722 b_fcus->nni->lk0 = lk0;
723 b_fcus->nni->lk1 = lk1;
724 b_fcus->nni->lk2 = lk2;
725
726 b_fcus->nni->score = lk0 - MAX(lk1,lk2);
727
728 if((b_fcus->nni->score < tree->mod->s_opt->min_diff_lk_local) &&
729 (b_fcus->nni->score > -tree->mod->s_opt->min_diff_lk_local))
730 {
731 b_fcus->nni->score = .0;
732 b_fcus->nni->lk1 = b_fcus->nni->lk2 = b_fcus->nni->lk0;
733 }
734
735
736 if((b_fcus->nni->lk1 > b_fcus->nni->lk0) && (b_fcus->nni->lk1 > b_fcus->nni->lk2))
737 {
738 if(b_fcus->nni->lk0 > b_fcus->nni->lk2) result = 1; //lk1 > lk0 > lk2
739 else result = 3; //lk1 > lk2 > lk0
740 }
741 else if((b_fcus->nni->lk2 > b_fcus->nni->lk0) && (b_fcus->nni->lk2 > b_fcus->nni->lk1))
742 {
743 if(b_fcus->nni->lk0 > b_fcus->nni->lk1) result = 2; //lk2 > lk0 > lk1
744 else result = 4; //lk2 > lk1 > lk0
745 }
746
747
748 Free_Scalar_Dbl(len_e1);
749 Free_Scalar_Dbl(len_e2);
750 Free_Scalar_Dbl(len_e3);
751 Free_Scalar_Dbl(len_e4);
752 Free_Scalar_Dbl(l_init);
753 Free_Scalar_Dbl(var_e1);
754 Free_Scalar_Dbl(var_e2);
755 Free_Scalar_Dbl(var_e3);
756 Free_Scalar_Dbl(var_e4);
757 Free_Scalar_Dbl(v_init);
758
759 return result;
760 }
761
762 //////////////////////////////////////////////////////////////
763 //////////////////////////////////////////////////////////////
764
765 /*
766 * Make one target swap, optimizing five branches.
767 * param tree : the tree to check
768 * param tested_t_edge : the swaping t_edge of the tree
769 * param swapToDo : 1 or 2, to select the first or the second swap to do
770 */
771
Make_Target_Swap(t_tree * tree,t_edge * b_fcus,int swaptodo)772 void Make_Target_Swap(t_tree *tree, t_edge *b_fcus, int swaptodo)
773 {
774 t_node *v1,*v2,*v3,*v4;
775 phydbl lktodo;
776 phydbl lk_init, lk_temp;
777 int i;
778
779 if(swaptodo < 0)
780 {
781 PhyML_Fprintf(stderr,"\n== Err in file %s at line %d\n\n",__FILE__,__LINE__);
782 Warn_And_Exit("");
783 }
784
785 //Initialization
786 lk_init = tree->c_lnL;
787
788 b_fcus->nni->score = .0;
789
790 lktodo = UNLIKELY;
791 v1 = v2 = v3 = v4 = NULL;
792
793 v1 = b_fcus->left->v[b_fcus->l_v1];
794 v2 = b_fcus->left->v[b_fcus->l_v2];
795 v3 = b_fcus->rght->v[b_fcus->r_v1];
796 v4 = b_fcus->rght->v[b_fcus->r_v2];
797
798 if(v1->num < v2->num)
799 {
800 PhyML_Fprintf(stderr,"\n== Err in file %s at line %d\n\n",__FILE__,__LINE__);
801 Warn_And_Exit("");
802 }
803 if(v3->num < v4->num)
804 {
805 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
806 Warn_And_Exit("");
807 }
808
809
810 /***********/
811 //make the selected swap
812 if(swaptodo==1)
813 {
814 Swap(v2,b_fcus->left,b_fcus->rght,v3,tree);
815 if(!Check_Topo_Constraints(tree,tree->io->cstr_tree)) Swap(v3,b_fcus->left,b_fcus->rght,v2,tree);
816 }
817 else
818 {
819 Swap(v2,b_fcus->left,b_fcus->rght,v4,tree);
820 if(!Check_Topo_Constraints(tree,tree->io->cstr_tree)) Swap(v4,b_fcus->left,b_fcus->rght,v2,tree);
821 }
822
823 MIXT_Set_Alias_Subpatt(YES,tree);
824 Set_Both_Sides(YES,tree);
825 lktodo = Update_Lk_At_Given_Edge(b_fcus,tree);
826
827 for(i=0;i<3;i++)
828 if(b_fcus->left->v[i] != b_fcus->rght)
829 Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
830
831 for(i=0;i<3;i++)
832 if(b_fcus->rght->v[i] != b_fcus->left)
833 Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
834
835 MIXT_Set_Alias_Subpatt(NO,tree);
836
837 //Optimize branch lengths and update likelihoods
838 lk_temp = UNLIKELY;
839 do
840 {
841 lktodo = lk_temp;
842
843 for(i=0;i<3;i++)
844 if(b_fcus->left->v[i] != b_fcus->rght)
845 {
846 Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
847 lk_temp = Br_Len_Opt(&(b_fcus->left->b[i]->l->v),b_fcus->left->b[i],tree);
848 }
849
850
851 Update_Partial_Lk(tree,b_fcus,b_fcus->left);
852 lk_temp = Br_Len_Opt(&(b_fcus->l->v),b_fcus,tree);
853
854 for(i=0;i<3;i++)
855 if(b_fcus->rght->v[i] != b_fcus->left)
856 {
857 Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
858 lk_temp = Br_Len_Opt(&(b_fcus->rght->b[i]->l->v),b_fcus->rght->b[i],tree);
859 }
860
861 Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
862
863 if(lk_temp < lktodo - tree->mod->s_opt->min_diff_lk_local)
864 {
865 PhyML_Fprintf(stderr,"\n== Edge %3d lk_temp = %f lktodo = %f\n",b_fcus->num,lk_temp,lktodo);
866 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
867 Warn_And_Exit("");
868 }
869 }
870 while(FABS(lk_temp-lktodo) > tree->mod->s_opt->min_diff_lk_global);
871 //until no significative improvement is detected
872
873
874 /* PhyML_Printf("\n.<< [%3d] l=%f lk_init=%f tree->c_lnL=%f score=%12f v1=%3d v2=%3d v3=%3d v4=%3d >>", */
875 /* b_fcus->num, */
876 /* b_fcus->l->v, */
877 /* lk_init, */
878 /* tree->c_lnL, */
879 /* lk_init-tree->c_lnL, */
880 /* v1->num,v2->num,v3->num,v4->num); */
881
882
883 if(tree->c_lnL < lk_init - tree->mod->s_opt->min_diff_lk_global)
884 {
885 PhyML_Fprintf(stderr,"\n== [%3d] v1=%d v2=%d v3=%d v4=%d",b_fcus->num,v1->num,v2->num,v3->num,v4->num);
886 PhyML_Fprintf(stderr,"\n== tree->c_lnL = %f lk_init = %f\n",tree->c_lnL,lk_init);
887 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
888 Warn_And_Exit("");
889 }
890 }
891
892 //////////////////////////////////////////////////////////////
893 //////////////////////////////////////////////////////////////
894
895 /**
896 * Convert an aLRT statistic to a non-parametric support
897 * param in: the statistic
898 */
899
Statistics_To_Probabilities(phydbl in)900 phydbl Statistics_To_Probabilities(phydbl in)
901 {
902 phydbl rough_value=0.0;
903 phydbl a=0.0;
904 phydbl b=0.0;
905 phydbl fa=0.0;
906 phydbl fb=0.0;
907
908 if(in>=0.000000393 && in<0.00000157)
909 {
910 a=0.000000393;
911 b=0.00000157;
912 fa=0.0005;
913 fb=0.001;
914 }
915 else if(in>=0.00000157 && in<0.0000393)
916 {
917 a=0.00000157;
918 b=0.0000393;
919 fa=0.001;
920 fb=0.005;
921 }
922 else if(in>=0.0000393 && in<0.000157)
923 {
924 a=0.0000393;
925 b=0.000157;
926 fa=0.005;
927 fb=0.01;
928 }
929 else if(in>=0.000157 && in<0.000982)
930 {
931 a=0.000157;
932 b=0.000982;
933 fa=0.01;
934 fb=0.025;
935 }
936 else if(in>0.000982 && in<0.00393)
937 {
938 a=0.000982;
939 b=0.00393;
940 fa=0.025;
941 fb=0.05;
942 }
943 else if(in>=0.00393 && in<0.0158)
944 {
945 a=0.00393;
946 b=0.0158;
947 fa=0.05;
948 fb=0.1;
949 }
950 else if(in>=0.0158 && in<0.0642)
951 {
952 a=0.0158;
953 b=0.0642;
954 fa=0.1;
955 fb=0.2;
956 }
957 else if(in>=0.0642 && in<0.148)
958 {
959 a=0.0642;
960 b=0.148;
961 fa=0.2;
962 fb=0.3;
963 }
964 else if(in>=0.148 && in<0.275)
965 {
966 a=0.148;
967 b=0.275;
968 fa=0.3;
969 fb=0.4;
970 }
971 else if(in>=0.275 && in<0.455)
972 {
973 a=0.275;
974 b=0.455;
975 fa=0.4;
976 fb=0.5;
977 }
978 else if(in>=0.455 && in<0.708)
979 {
980 a=0.455;
981 b=0.708;
982 fa=0.5;
983 fb=0.6;
984 }
985 else if(in>=0.708 && in<1.074)
986 {
987 a=0.708;
988 b=1.074;
989 fa=0.6;
990 fb=0.7;
991 }
992 else if(in>=1.074 && in<1.642)
993 {
994 a=1.074;
995 b=1.642;
996 fa=0.7;
997 fb=0.8;
998 }
999 else if(in>=1.642 && in<2.706)
1000 {
1001 a=1.642;
1002 b=2.706;
1003 fa=0.8;
1004 fb=0.9;
1005 }
1006 else if(in>=2.706 && in<3.841)
1007 {
1008 a=2.706;
1009 b=3.841;
1010 fa=0.9;
1011 fb=0.95;
1012 }
1013 else if(in>=3.841 && in<5.024)
1014 {
1015 a=3.841;
1016 b=5.024;
1017 fa=0.95;
1018 fb=0.975;
1019 }
1020 else if(in>=5.024 && in<6.635)
1021 {
1022 a=5.024;
1023 b=6.635;
1024 fa=0.975;
1025 fb=0.99;
1026 }
1027 else if(in>=6.635 && in<7.879)
1028 {
1029 a=6.635;
1030 b=7.879;
1031 fa=0.99;
1032 fb=0.995;
1033 }
1034 else if(in>=7.879 && in<10.828)
1035 {
1036 a=7.879;
1037 b=10.828;
1038 fa=0.995;
1039 fb=0.999;
1040 }
1041 else if(in>=10.828 && in<12.116)
1042 {
1043 a=10.828;
1044 b=12.116;
1045 fa=0.999;
1046 fb=0.9995;
1047 }
1048 if (in>=12.116)
1049 {
1050 rough_value=0.9999;
1051 }
1052 else if(in<0.000000393)
1053 {
1054 rough_value=0.0001;
1055 }
1056 else
1057 {
1058 rough_value=(b-in)/(b-a)*fa + (in - a)/(b-a)*fb;
1059 }
1060 rough_value=rough_value+(1.0-rough_value)/2.0;
1061 rough_value=rough_value*rough_value*rough_value;
1062 return rough_value;
1063 }
1064
1065 //////////////////////////////////////////////////////////////
1066 //////////////////////////////////////////////////////////////
1067
1068 /**
1069 * deprecated
1070 * Compute a RELL support, using the latest tested branch
1071 * param tree: the tested tree
1072 */
Statistics_to_RELL(t_tree * tree)1073 phydbl Statistics_to_RELL(t_tree *tree)
1074 {
1075 int i;
1076 int occurence=10000;
1077 phydbl nb=0.0;
1078 phydbl res,*pi;
1079 int site;
1080 phydbl lk0=0.0;
1081 phydbl lk1=0.0;
1082 phydbl lk2=0.0;
1083 int position = -1;
1084 t_tree *buff_tree;
1085 int* samples;
1086
1087 /*! 1000 times */
1088 for(i=0;i<occurence;i++)
1089 {
1090 lk0=0.0;
1091 lk1=0.0;
1092 lk2=0.0;
1093
1094 /*! Shuffle the data and increment the support, if needed */
1095 buff_tree = tree;
1096 do
1097 {
1098 pi = (phydbl *)mCalloc(tree->data->crunch_len,sizeof(phydbl));
1099 for(site=0;site<tree->n_pattern;site++) pi[site] = tree->data->wght[site] / (phydbl)tree->data->init_len;
1100
1101 samples = Sample_n_i_With_Proba_pi(pi,tree->n_pattern,tree->data->init_len);
1102 For(site, tree->data->init_len)
1103 {
1104 position = samples[site];
1105 lk0+=tree->log_lks_aLRT[0][position];
1106 lk1+=tree->log_lks_aLRT[1][position];
1107 lk2+=tree->log_lks_aLRT[2][position];
1108 }
1109 if (lk0>=lk1 && lk0>=lk2) nb++;
1110 tree = tree->next_mixt;
1111 Free(samples);
1112 Free(pi);
1113 }
1114 while(tree);
1115 tree = buff_tree;
1116 }
1117
1118 res= nb/(phydbl)occurence;
1119
1120 return res;
1121 }
1122 //////////////////////////////////////////////////////////////
1123 //////////////////////////////////////////////////////////////
1124
1125 /*!
1126 Compute a SH-like support, using the latest tested branch
1127 param tree: the tested tree
1128 */
Statistics_To_SH(t_tree * tree)1129 phydbl Statistics_To_SH(t_tree *tree)
1130 {
1131 int i;
1132 int occurence=10000;
1133 phydbl nb=0.0;
1134 phydbl res;
1135 int site;
1136 phydbl lk0=0.0;
1137 phydbl lk1=0.0;
1138 phydbl lk2=0.0;
1139 phydbl c0=0.0;
1140 phydbl c1=0.0;
1141 phydbl c2=0.0;
1142 int position=-1;
1143 phydbl delta_local=-1.;
1144 phydbl delta=0.0;
1145 t_tree *buff_tree;
1146 phydbl *pi;
1147 int* samples;
1148
1149 /*! Compute the total log-lk of each NNI position */
1150 buff_tree = tree;
1151 do
1152 {
1153 For(site, tree->n_pattern)
1154 {
1155 c0+=tree->log_lks_aLRT[0][site] * tree->data->wght[site];
1156 c1+=tree->log_lks_aLRT[1][site] * tree->data->wght[site];
1157 c2+=tree->log_lks_aLRT[2][site] * tree->data->wght[site];
1158 }
1159 tree = tree->next_mixt;
1160 }
1161 while(tree);
1162 tree = buff_tree;
1163
1164
1165 if (c0>=c1 && c0>=c2)
1166 {
1167 if (c1>=c2)
1168 {
1169 delta=c0-c1;
1170 }
1171 else
1172 {
1173 delta=c0-c2;
1174 }
1175 }
1176 else if(c1>=c0 && c1>=c2)
1177 {
1178 if (c0>=c2)
1179 {
1180 delta=c1-c0;
1181 }
1182 else
1183 {
1184 delta=c1-c2;
1185 }
1186 }
1187 else
1188 {
1189 if (c1>=c0)
1190 {
1191 delta=c2-c1;
1192 }
1193 else
1194 {
1195 delta=c2-c0;
1196 }
1197 }
1198
1199 /*! 1000 times */
1200 for(i=0;i<occurence;i++)
1201 {
1202 lk0=0.0;
1203 lk1=0.0;
1204 lk2=0.0;
1205
1206 buff_tree = tree;
1207 do
1208 {
1209 pi = (phydbl *)mCalloc(tree->data->crunch_len,sizeof(phydbl));
1210 for(site=0;site<tree->n_pattern;site++) pi[site] = tree->data->wght[site] / (phydbl)tree->data->init_len;
1211
1212 samples = Sample_n_i_With_Proba_pi(pi,tree->n_pattern,tree->data->init_len);
1213 /*! Shuffle the data */
1214 For(site, tree->data->init_len)
1215 {
1216 position = samples[site];
1217 lk0+=tree->log_lks_aLRT[0][position];
1218 lk1+=tree->log_lks_aLRT[1][position];
1219 lk2+=tree->log_lks_aLRT[2][position];
1220 }
1221
1222 tree = tree->next_mixt;
1223 Free(samples);
1224 Free(pi);
1225 }
1226 while(tree);
1227 tree = buff_tree;
1228
1229 /*! return to null hypothesis */
1230 lk0=lk0-c0;
1231 lk1=lk1-c1;
1232 lk2=lk2-c2;
1233
1234 /*! compute results and increment if needed */
1235 delta_local=0.0;
1236 if (lk0>=lk1 && lk0>=lk2)
1237 {
1238 if (lk1>=lk2)
1239 {
1240 delta_local=lk0-lk1;
1241 }
1242 else
1243 {
1244 delta_local=lk0-lk2;
1245 }
1246 }
1247 else if(lk1>=lk0 && lk1>=lk2)
1248 {
1249 if (lk0>=lk2)
1250 {
1251 delta_local=lk1-lk0;
1252 }
1253 else
1254 {
1255 delta_local=lk1-lk2;
1256 }
1257 }
1258 else
1259 {
1260 if (lk1>=lk0)
1261 {
1262 delta_local=lk2-lk1;
1263 }
1264 else
1265 {
1266 delta_local=lk2-lk0;
1267 }
1268 }
1269
1270 if (delta>(delta_local+0.1))
1271 {
1272 nb++;
1273 }
1274 }
1275
1276 res= nb/occurence;
1277
1278 return res;
1279 }
1280
1281 //////////////////////////////////////////////////////////////
1282 //////////////////////////////////////////////////////////////
1283
1284 /**
1285 * deprecated
1286 * Compute one side likelihood
1287 * param b_fcus : concerned edge
1288 * param tree : b_fcus tree
1289 * param exclude : side to exclude for computation
1290 */
Update_Lk_At_Given_Edge_Excluding(t_edge * b_fcus,t_tree * tree,t_node * exclude)1291 phydbl Update_Lk_At_Given_Edge_Excluding(t_edge *b_fcus, t_tree *tree, t_node *exclude)
1292 {
1293 if((!b_fcus->left->tax) && (exclude == NULL || exclude != b_fcus->left))
1294 Update_Partial_Lk(tree,b_fcus,b_fcus->left);
1295 if((!b_fcus->rght->tax) && (exclude == NULL || exclude != b_fcus->rght))
1296 Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
1297
1298 tree->c_lnL = Lk(b_fcus,tree);
1299
1300 return tree->c_lnL;
1301 }
1302
1303
1304 //////////////////////////////////////////////////////////////
1305 //////////////////////////////////////////////////////////////
1306
1307 //////////////////////////////////////////////////////////////
1308 //////////////////////////////////////////////////////////////
1309
1310 //////////////////////////////////////////////////////////////
1311 //////////////////////////////////////////////////////////////
1312
1313 //////////////////////////////////////////////////////////////
1314 //////////////////////////////////////////////////////////////
1315
1316 //////////////////////////////////////////////////////////////
1317 //////////////////////////////////////////////////////////////
1318
1319 //////////////////////////////////////////////////////////////
1320 //////////////////////////////////////////////////////////////
1321
1322