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 "simu.h"
14 #ifdef BEAGLE
15 #include "beagle_utils.h"
16 #endif
17
18
19 //////////////////////////////////////////////////////////////
20 //////////////////////////////////////////////////////////////
21
Simu_Loop(t_tree * tree)22 void Simu_Loop(t_tree *tree)
23 {
24 return;
25 }
26
27 //////////////////////////////////////////////////////////////
28 //////////////////////////////////////////////////////////////
29
Simu(t_tree * tree,int n_step_max,phydbl delta_lnL,phydbl init_T,phydbl delta_T,int min_n_edges_traversed)30 int Simu(t_tree *tree, int n_step_max, phydbl delta_lnL, phydbl init_T, phydbl delta_T, int min_n_edges_traversed)
31 {
32 phydbl old_loglk,delta;
33 unsigned int n_round;
34 time_t t_cur;
35
36 tree->c_lnL = UNLIKELY;
37 delta = UNLIKELY;
38 old_loglk = UNLIKELY;
39 n_round = 0;
40 tree->annealing_temp = init_T;
41 do
42 {
43 for(int i=0;i<2*tree->n_otu-3;++i) tree->a_edges[i]->l->v *= Rgamma((phydbl)(0.2*n_round+1),(phydbl)(1./(0.2*n_round+1)));
44 old_loglk = tree->c_lnL;
45 Set_Both_Sides(NO,tree);
46 tree->tip_root = Rand_Int(0,tree->n_otu-1);
47 Lk(NULL,tree);
48 tree->n_edges_traversed = 0;
49 tree->fully_nni_opt = YES;
50 NNI_Traversal(tree->a_nodes[tree->tip_root],
51 tree->a_nodes[tree->tip_root]->v[0],
52 NULL,
53 tree->a_nodes[tree->tip_root]->b[0],
54 YES,
55 tree);
56 delta = tree->c_lnL - old_loglk;
57 tree->annealing_temp -= delta_T;
58 if(tree->annealing_temp < 0.0) tree->annealing_temp = 0.0;
59 n_round++;
60 time(&t_cur);
61
62 PhyML_Printf("\n. %5ds lnL: %12G T: %12G %4d/%4d",
63 (int)(t_cur-tree->t_beg),
64 tree->c_lnL,
65 tree->annealing_temp,
66 tree->n_edges_traversed,
67 tree->n_otu);
68
69 if((n_round >= n_step_max || tree->fully_nni_opt == YES) && Are_Equal(tree->annealing_temp,0.0,1.E-3) && delta < delta_lnL) break;
70 }
71 while(1);
72
73 return 1;
74 }
75
76 //////////////////////////////////////////////////////////////
77 //////////////////////////////////////////////////////////////
78
79
Simu_Pars(t_tree * tree,int n_step_max)80 void Simu_Pars(t_tree *tree, int n_step_max)
81 {
82 phydbl old_pars,n_iter,lambda;
83 int i,n_neg,n_tested,n_without_swap,n_tot_swap,step;
84 t_edge **sorted_b,**tested_b;
85 int each;
86
87 sorted_b = (t_edge **)mCalloc(tree->n_otu-3,sizeof(t_edge *));
88 tested_b = (t_edge **)mCalloc(tree->n_otu-3,sizeof(t_edge *));
89
90 old_pars = 0;
91 tree->c_pars = 0;
92 n_iter = 1.0;
93 n_tested = 0;
94 n_without_swap = 0;
95 step = 0;
96 each = 4;
97 lambda = 0.5;
98 n_tot_swap = 0;
99
100 Update_Dirs(tree);
101
102 if((tree->verbose > VL0) && (tree->io->quiet == NO)) PhyML_Printf("\n\n. Starting simultaneous NNI moves (parsimony criterion)...\n");
103
104 do
105 {
106 ++step;
107
108 if(step > n_step_max) break;
109
110 each--;
111
112 Set_Both_Sides(YES,tree);
113 Pars(NULL,tree);
114
115 if((tree->verbose > VL0) && (tree->io->quiet == NO))
116 {
117 Print_Pars(tree);
118 if(step > 1) (n_tested > 1)?(printf("[%4d NNIs]",n_tested)):(printf("[%4d NNI ]",n_tested));
119 }
120
121
122 if(FABS(old_pars - tree->c_pars) < SMALL) break;
123
124 if((tree->c_pars > old_pars) && (step > 1))
125 {
126 if(tree->verbose > VL0 && tree->io->quiet == NO)
127 PhyML_Printf("\n\n. Moving backward (topology) \n");
128 if(!Mov_Backward_Topo_Pars(tree,old_pars,tested_b,n_tested))
129 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
130 if(!tree->n_swap) n_neg = 0;
131
132 Set_Both_Sides(YES,tree);
133 Pars(NULL,tree);
134 }
135 else
136 {
137 old_pars = tree->c_pars;
138
139 n_neg = 0;
140 For(i,2*tree->n_otu-3)
141 if((!tree->a_edges[i]->left->tax) &&
142 (!tree->a_edges[i]->rght->tax))
143 NNI_Pars(tree,tree->a_edges[i],NO);
144
145 Select_Edges_To_Swap(tree,sorted_b,&n_neg);
146 Sort_Edges_NNI_Score(tree,sorted_b,n_neg);
147
148 n_tested = 0;
149 For(i,(int)CEIL((phydbl)n_neg*(lambda)))
150 tested_b[n_tested++] = sorted_b[i];
151
152 Make_N_Swap(tree,tested_b,0,n_tested);
153
154 n_tot_swap += n_tested;
155
156 if(n_tested > 0) n_without_swap = 0;
157 else n_without_swap++;
158 }
159 n_iter+=1.0;
160 }
161 while(1);
162
163 Free(sorted_b);
164 Free(tested_b);
165 }
166
167 //////////////////////////////////////////////////////////////
168 //////////////////////////////////////////////////////////////
169
Select_Edges_To_Swap(t_tree * tree,t_edge ** sorted_b,int * n_neg)170 void Select_Edges_To_Swap(t_tree *tree, t_edge **sorted_b, int *n_neg)
171 {
172 int i;
173 t_edge *b;
174 /* phydbl best_score; */
175
176 *n_neg = 0;
177
178 For(i,2*tree->n_otu-3)
179 {
180 b = tree->a_edges[i];
181 /* best_score = b->nni->score; */
182
183 if((!b->left->tax) && (!b->rght->tax) && (b->nni->score < -tree->mod->s_opt->min_diff_lk_move))
184 {
185 /* // Evaluate NNIs on edges at distance 1 */
186 /* Check_NNI_Scores_Around(b->left,b->rght,b,&best_score,tree); */
187 /* Check_NNI_Scores_Around(b->rght,b->left,b,&best_score,tree); */
188
189 /* // Evaluate NNIs on edges at distance 2 */
190 /* Check_NNI_Scores_Around(b->left,b->left->v[b->l_v1],b,&best_score,tree); */
191 /* Check_NNI_Scores_Around(b->left,b->left->v[b->l_v2],b,&best_score,tree); */
192
193 /* // Evaluate NNIs on edges at distance 2 */
194 /* Check_NNI_Scores_Around(b->rght,b->rght->v[b->r_v1],b,&best_score,tree); */
195 /* Check_NNI_Scores_Around(b->rght,b->rght->v[b->r_v2],b,&best_score,tree); */
196
197 /* if(best_score < b->nni->score) continue; */
198
199 sorted_b[*n_neg] = b;
200 (*n_neg)++;
201 }
202
203
204 /* if((!b->left->tax) && (!b->rght->tax) && (b->nni->score < -tree->mod->s_opt->min_diff_lk_move)) */
205 /* { */
206 /* sorted_b[*n_neg] = b; */
207 /* (*n_neg)++; */
208 /* } */
209
210 }
211 }
212
213 //////////////////////////////////////////////////////////////
214 //////////////////////////////////////////////////////////////
215
216
Update_Bl(t_tree * tree,phydbl fact)217 void Update_Bl(t_tree *tree, phydbl fact)
218 {
219 int i;
220 scalar_dbl *l,*l_old,*l0;
221
222 For(i,2*tree->n_otu-3)
223 {
224 l = tree->a_edges[i]->l;
225 l_old = tree->a_edges[i]->l_old;
226 l0 = tree->a_edges[i]->nni->l0;
227
228 do
229 {
230 l->v = l_old->v + (l0->v - l_old->v)*fact;
231 l = l->next;
232 l_old = l_old->next;
233 l0 = l0->next;
234 }
235 while(l);
236 }
237 }
238
239 //////////////////////////////////////////////////////////////
240 //////////////////////////////////////////////////////////////
241
Make_N_Swap(t_tree * tree,t_edge ** b,int beg,int end)242 void Make_N_Swap(t_tree *tree,t_edge **b, int beg, int end)
243 {
244 int i;
245 /* t_edge *orig; */
246 t_node *n1,*n2,*n3,*n4;
247
248 n1 = n2 = n3 = n4 = NULL;
249
250 tree->n_swap = 0;
251 for(i=beg;i<end;i++)
252 {
253 n1 = n2 = n3 = n4 = NULL;
254
255 n1 = b[i]->nni->swap_node_v1;
256 n2 = b[i]->nni->swap_node_v2;
257 n3 = b[i]->nni->swap_node_v3;
258 n4 = b[i]->nni->swap_node_v4;
259
260 /* if(b[i]->nni->best_conf == 1) */
261 /* { */
262 /* n1 = b[i]->left->v[b[i]->l_v2]; */
263 /* n2 = b[i]->left; */
264 /* n3 = b[i]->rght; */
265 /* n4 = b[i]->rght->v[b[i]->r_v1]; */
266 /* } */
267 /* else if(b[i]->nni->best_conf == 2) */
268 /* { */
269 /* n1 = b[i]->left->v[b[i]->l_v2]; */
270 /* n2 = b[i]->left; */
271 /* n3 = b[i]->rght; */
272 /* n4 = b[i]->rght->v[b[i]->r_v2]; */
273 /* } */
274
275 if(b[i]->nni->best_conf == 1)
276 {
277 if(n1 != b[i]->left->v[b[i]->l_v2] ||
278 /* n2 != b[i]->left || */
279 /* n3 != b[i]->rght || */
280 n4 != b[i]->rght->v[b[i]->r_v1]) continue;
281 }
282 else if(b[i]->nni->best_conf == 2)
283 {
284 if(n1 != b[i]->left->v[b[i]->l_v2] ||
285 /* n2 != b[i]->left || */
286 /* n3 != b[i]->rght || */
287 n4 != b[i]->rght->v[b[i]->r_v2]) continue;
288 }
289
290 Swap(n1,n2,n3,n4,tree);
291
292 if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
293 {
294 /* Undo this swap as it violates one of the topological constraints
295 defined in the input constraint tree
296 */
297 Swap(n4,n2,n3,n1,tree);
298 }
299
300 if(tree->n_root)
301 {
302 tree->n_root->v[2] = tree->e_root->left;
303 tree->n_root->v[1] = tree->e_root->rght;
304 }
305
306 Copy_Scalar_Dbl(b[i]->nni->best_l,b[i]->l);
307 Copy_Scalar_Dbl(b[i]->nni->best_v,b[i]->l_var);
308
309 /* orig = b[i]; */
310 /* do */
311 /* { */
312 /* b[i]->l->v = b[i]->nni->best_l; */
313 /* b[i] = b[i]->next; */
314 /* } */
315 /* while(b[i]); */
316 /* b[i] = orig; */
317
318 tree->n_swap++;
319 }
320
321
322 /* PhyML_Printf("\n. End Actually performing swaps\n"); */
323
324 }
325
326 //////////////////////////////////////////////////////////////
327 //////////////////////////////////////////////////////////////
328
Make_Best_Swap(t_tree * tree)329 int Make_Best_Swap(t_tree *tree)
330 {
331 int i,j,return_value;
332 t_edge *b,**sorted_b;
333 /* t_edge *orig; */
334 t_node *n1,*n2,*n3,*n4;
335
336
337 sorted_b = (t_edge **)mCalloc(tree->n_otu-3,sizeof(t_edge *));
338
339 j=0;
340 For(i,2*tree->n_otu-3)
341 if((!tree->a_edges[i]->left->tax) &&
342 (!tree->a_edges[i]->rght->tax))
343 sorted_b[j++] = tree->a_edges[i];
344
345 Sort_Edges_NNI_Score(tree,sorted_b,tree->n_otu-3);
346
347 if(sorted_b[0]->nni->score < -0.0)
348 {
349 b = sorted_b[0];
350 return_value = 1;
351
352 n1 = n2 = n3 = n4 = NULL;
353
354 /* n1 = b->nni->swap_node_v1; */
355 /* n2 = b->nni->swap_node_v2; */
356 /* n3 = b->nni->swap_node_v3; */
357 /* n4 = b->nni->swap_node_v4; */
358
359 if(b->nni->best_conf == 1)
360 {
361 n1 = b->left->v[b->l_v2];
362 n2 = b->left;
363 n3 = b->rght;
364 n4 = b->rght->v[b->r_v1];
365 }
366 else if(b->nni->best_conf == 2)
367 {
368 n1 = b->left->v[b->l_v2];
369 n2 = b->left;
370 n3 = b->rght;
371 n4 = b->rght->v[b->r_v2];
372 }
373
374 Swap(n1,n2,n3,n4,tree);
375
376 if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
377 {
378 /* Undo this swap as it violates one of the topological constraints
379 defined in the input constraint tree
380 */
381 Swap(n4,n2,n3,n1,tree);
382 }
383
384 /* b->l->v = b->nni->best_l; */
385
386 Copy_Scalar_Dbl(b->nni->best_l,b->l);
387 Copy_Scalar_Dbl(b->nni->best_v,b->l_var);
388
389 /* orig = b; */
390 /* do */
391 /* { */
392 /* b->l->v = b->nni->best_l; */
393 /* b = b->next; */
394 /* } */
395 /* while(b); */
396 /* b = orig; */
397 }
398 else return_value = 0;
399
400 Free(sorted_b);
401
402 return return_value;
403 }
404
405 //////////////////////////////////////////////////////////////
406 //////////////////////////////////////////////////////////////
407
Mov_Backward_Topo_Bl(t_tree * tree,phydbl lk_old,t_edge ** tested_b,int n_tested)408 int Mov_Backward_Topo_Bl(t_tree *tree, phydbl lk_old, t_edge **tested_b, int n_tested)
409 {
410 scalar_dbl **l_init,**v_init;
411 int i,step,beg,end;
412 t_edge *b,*orig;
413
414 l_init = (scalar_dbl **)mCalloc(2*tree->n_otu-3,sizeof(scalar_dbl *));
415 v_init = (scalar_dbl **)mCalloc(2*tree->n_otu-3,sizeof(scalar_dbl *));
416
417 For(i,2*tree->n_otu-3)
418 {
419 l_init[i] = Duplicate_Scalar_Dbl(tree->a_edges[i]->l);
420 v_init[i] = Duplicate_Scalar_Dbl(tree->a_edges[i]->l_var);
421 }
422
423 step = 2;
424 do
425 {
426 For(i,2*tree->n_otu-3)
427 {
428 b = tree->a_edges[i];
429
430 orig = b;
431 do
432 {
433 b->l->v = b->l_old->v + (1./step) * (l_init[i]->v - b->l_old->v);
434 if(b->next == NULL) break;
435 b = b->next;
436 l_init[i] = l_init[i]->next;
437 }
438 while(b);
439 b = orig;
440 }
441
442 beg = (int)FLOOR((phydbl)n_tested/(step-1));
443 end = 0;
444 Unswap_N_Branch(tree,tested_b,beg,end);
445 beg = 0;
446 end = (int)FLOOR((phydbl)n_tested/step);
447 Swap_N_Branch(tree,tested_b,beg,end);
448
449 if(!end) tree->n_swap = 0;
450
451 Set_Both_Sides(NO,tree);
452 Lk(NULL,tree);
453
454 step++;
455
456 }while((tree->c_lnL < lk_old) && (step < 1000));
457
458
459 if(step == 1000)
460 {
461 if(tree->n_swap) Exit("\n== Err. in Mov_Backward_Topo_Bl (n_swap > 0)\n");
462
463 Restore_Br_Len(tree);
464
465 Set_Both_Sides(NO,tree);
466 Lk(NULL,tree);
467 }
468
469 For(i,2*tree->n_otu-3)
470 {
471 while(l_init[i]->prev) l_init[i] = l_init[i]->prev;
472 Free_Scalar_Dbl(l_init[i]);
473 }
474 Free(l_init);
475
476 For(i,2*tree->n_otu-3)
477 {
478 while(v_init[i]->prev) v_init[i] = v_init[i]->prev;
479 Free_Scalar_Dbl(v_init[i]);
480 }
481 Free(v_init);
482
483 tree->n_swap = 0;
484 For(i,2*tree->n_otu-3)
485 {
486 if(tree->a_edges[i]->nni->score < 0.0) tree->n_swap++;
487 tree->a_edges[i]->nni->score = +1.0;
488 }
489
490
491 if(tree->c_lnL > lk_old) return 1;
492 else if((tree->c_lnL > lk_old-tree->mod->s_opt->min_diff_lk_local) &&
493 (tree->c_lnL < lk_old+tree->mod->s_opt->min_diff_lk_local)) return -1;
494 else return 0;
495 }
496
497 //////////////////////////////////////////////////////////////
498 //////////////////////////////////////////////////////////////
499
Mov_Backward_Topo_Pars(t_tree * tree,int pars_old,t_edge ** tested_b,int n_tested)500 int Mov_Backward_Topo_Pars(t_tree *tree, int pars_old, t_edge **tested_b, int n_tested)
501 {
502 int i,step,beg,end;
503
504 step = 2;
505 do
506 {
507 beg = (int)FLOOR((phydbl)n_tested/(step-1));
508 end = 0;
509 Unswap_N_Branch(tree,tested_b,beg,end);
510 beg = 0;
511 end = (int)FLOOR((phydbl)n_tested/step);
512 Swap_N_Branch(tree,tested_b,beg,end);
513
514 if(!end) tree->n_swap = 0;
515
516 Set_Both_Sides(NO,tree);
517 Pars(NULL,tree);
518
519 step++;
520
521 }
522 while((tree->c_pars > pars_old) && (step < 1000));
523
524
525 if(step == 1000)
526 {
527 if(tree->n_swap) Exit("\n. Err. in Mov_Backward_Topo_Bl (n_swap > 0)\n");
528
529 Set_Both_Sides(NO,tree);
530 Pars(NULL,tree);
531 }
532
533 tree->n_swap = 0;
534 For(i,2*tree->n_otu-3)
535 {
536 if(tree->a_edges[i]->nni->score < 0.0) tree->n_swap++;
537 tree->a_edges[i]->nni->score = +1.0;
538 }
539
540
541 if(tree->c_pars < pars_old) return 1;
542 else if(tree->c_pars == pars_old) return -1;
543 else return 0;
544 }
545
546 //////////////////////////////////////////////////////////////
547 //////////////////////////////////////////////////////////////
548
549
Unswap_N_Branch(t_tree * tree,t_edge ** b,int beg,int end)550 void Unswap_N_Branch(t_tree *tree, t_edge **b, int beg, int end)
551 {
552 int i;
553 /* t_edge *orig; */
554 t_node *n1,*n2,*n3,*n4;
555
556 n1 = n2 = n3 = n4 = NULL;
557
558 if(end>beg)
559 {
560 for(i=beg;i<end;i++)
561 {
562 n1 = n2 = n3 = n4 = NULL;
563
564 n1 = b[i]->nni->swap_node_v4;
565 n2 = b[i]->nni->swap_node_v2;
566 n3 = b[i]->nni->swap_node_v3;
567 n4 = b[i]->nni->swap_node_v1;
568
569 /* if(b[i]->nni->best_conf == 1) */
570 /* { */
571 /* n1 = b[i]->left->v[b[i]->l_v2]; */
572 /* n2 = b[i]->left; */
573 /* n3 = b[i]->rght; */
574 /* n4 = b[i]->rght->v[b[i]->r_v1]; */
575 /* } */
576 /* else if(b[i]->nni->best_conf == 2) */
577 /* { */
578 /* n1 = b[i]->left->v[b[i]->l_v2]; */
579 /* n2 = b[i]->left; */
580 /* n3 = b[i]->rght; */
581 /* n4 = b[i]->rght->v[b[i]->r_v2]; */
582 /* } */
583
584 Swap(n1,n2,n3,n4,tree);
585
586 if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
587 {
588 /* Undo this swap as it violates one of the topological constraints
589 defined in the input constraint tree
590 */
591 Swap(n4,n2,n3,n1,tree);
592 }
593
594
595 /* (b[i]->nni->best_conf == 1)? */
596 /* (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1],tree)): */
597 /* (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2],tree)); */
598
599 /* b[i]->l->v = b[i]->l_old->v; */
600
601 Copy_Scalar_Dbl(b[i]->l_old,b[i]->l);
602 Copy_Scalar_Dbl(b[i]->l_var_old,b[i]->l_var);
603
604 /* orig = b[i]; */
605 /* do */
606 /* { */
607 /* b[i]->l->v = b[i]->l_old->v; */
608 /* b[i] = b[i]->next; */
609 /* } */
610 /* while(b[i]); */
611 /* b[i] = orig; */
612
613 }
614 }
615 else
616 {
617 for(i=beg-1;i>=end;i--)
618 {
619 n1 = n2 = n3 = n4 = NULL;
620
621 n1 = b[i]->nni->swap_node_v4;
622 n2 = b[i]->nni->swap_node_v2;
623 n3 = b[i]->nni->swap_node_v3;
624 n4 = b[i]->nni->swap_node_v1;
625
626 /* if(b[i]->nni->best_conf == 1) */
627 /* { */
628 /* n1 = b[i]->left->v[b[i]->l_v2]; */
629 /* n2 = b[i]->left; */
630 /* n3 = b[i]->rght; */
631 /* n4 = b[i]->rght->v[b[i]->r_v1]; */
632 /* } */
633 /* else if(b[i]->nni->best_conf == 2) */
634 /* { */
635 /* n1 = b[i]->left->v[b[i]->l_v2]; */
636 /* n2 = b[i]->left; */
637 /* n3 = b[i]->rght; */
638 /* n4 = b[i]->rght->v[b[i]->r_v2]; */
639 /* } */
640
641 Swap(n1,n2,n3,n4,tree);
642
643 if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
644 {
645 /* Undo this swap as it violates one of the topological constraints
646 defined in the input constraint tree
647 */
648 Swap(n4,n2,n3,n1,tree);
649 }
650
651
652 /* b[i]->l->v = b[i]->l_old->v; */
653 Copy_Scalar_Dbl(b[i]->l_old,b[i]->l);
654 Copy_Scalar_Dbl(b[i]->l_var_old,b[i]->l_var);
655
656
657 /* orig = b[i]; */
658 /* do */
659 /* { */
660 /* b[i]->l->v = b[i]->l_old->v; */
661 /* b[i] = b[i]->next; */
662 /* } */
663 /* while(b[i]); */
664 /* b[i] = orig; */
665
666 }
667 }
668 }
669
670 //////////////////////////////////////////////////////////////
671 //////////////////////////////////////////////////////////////
672
Swap_N_Branch(t_tree * tree,t_edge ** b,int beg,int end)673 void Swap_N_Branch(t_tree *tree,t_edge **b, int beg, int end)
674 {
675 int i;
676 /* t_edge *orig; */
677 t_node *n1,*n2,*n3,*n4;
678
679 n1 = n2 = n3 = n4 = NULL;
680
681 if(end>beg)
682 {
683 for(i=beg;i<end;i++)
684 {
685 n1 = n2 = n3 = n4 = NULL;
686
687 n1 = b[i]->nni->swap_node_v1;
688 n2 = b[i]->nni->swap_node_v2;
689 n3 = b[i]->nni->swap_node_v3;
690 n4 = b[i]->nni->swap_node_v4;
691
692
693 /* if(b[i]->nni->best_conf == 1) */
694 /* { */
695 /* n1 = b[i]->left->v[b[i]->l_v2]; */
696 /* n2 = b[i]->left; */
697 /* n3 = b[i]->rght; */
698 /* n4 = b[i]->rght->v[b[i]->r_v1]; */
699 /* } */
700 /* else if(b[i]->nni->best_conf == 2) */
701 /* { */
702 /* n1 = b[i]->left->v[b[i]->l_v2]; */
703 /* n2 = b[i]->left; */
704 /* n3 = b[i]->rght; */
705 /* n4 = b[i]->rght->v[b[i]->r_v2]; */
706 /* } */
707
708 Swap(n1,n2,n3,n4,tree);
709
710 if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
711 {
712 /* Undo this swap as it violates one of the topological constraints
713 defined in the input constraint tree
714 */
715 Swap(n4,n2,n3,n1,tree);
716 }
717
718 /* b[i]->l->v = b[i]->nni->best_l; */
719
720 Copy_Scalar_Dbl(b[i]->nni->best_l,b[i]->l);
721 Copy_Scalar_Dbl(b[i]->nni->best_v,b[i]->l_var);
722
723 /* orig = b[i]; */
724 /* do */
725 /* { */
726 /* b[i]->l->v = b[i]->nni->best_l; */
727 /* b[i] = b[i]->next; */
728 /* } */
729 /* while(b[i]); */
730 /* b[i] = orig; */
731 }
732 }
733 else
734 {
735 for(i=beg-1;i>=end;i--)
736 {
737
738 n1 = n2 = n3 = n4 = NULL;
739
740 n1 = b[i]->nni->swap_node_v1;
741 n2 = b[i]->nni->swap_node_v2;
742 n3 = b[i]->nni->swap_node_v3;
743 n4 = b[i]->nni->swap_node_v4;
744
745
746 /* if(b[i]->nni->best_conf == 1) */
747 /* { */
748 /* n1 = b[i]->left->v[b[i]->l_v2]; */
749 /* n2 = b[i]->left; */
750 /* n3 = b[i]->rght; */
751 /* n4 = b[i]->rght->v[b[i]->r_v1]; */
752 /* } */
753 /* else if(b[i]->nni->best_conf == 2) */
754 /* { */
755 /* n1 = b[i]->left->v[b[i]->l_v2]; */
756 /* n2 = b[i]->left; */
757 /* n3 = b[i]->rght; */
758 /* n4 = b[i]->rght->v[b[i]->r_v2]; */
759 /* } */
760
761 Swap(n1,n2,n3,n4,tree);
762
763 if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
764 {
765 /* Undo this swap as it violates one of the topological constraints
766 defined in the input constraint tree
767 */
768 Swap(n4,n2,n3,n1,tree);
769 }
770
771 /* b[i]->l->v = b[i]->nni->best_l; */
772
773 Copy_Scalar_Dbl(b[i]->nni->best_l,b[i]->l);
774 Copy_Scalar_Dbl(b[i]->nni->best_v,b[i]->l_var);
775
776 /* orig = b[i]; */
777 /* do */
778 /* { */
779 /* b[i]->l->v = b[i]->nni->best_l; */
780 /* b[i] = b[i]->next; */
781 /* } */
782 /* while(b[i]); */
783 /* b[i] = orig; */
784 }
785 }
786 }
787
788 //////////////////////////////////////////////////////////////
789 //////////////////////////////////////////////////////////////
790
Check_NNI_Scores_Around(t_node * a,t_node * d,t_edge * b,phydbl * best_score,t_tree * tree)791 void Check_NNI_Scores_Around(t_node *a, t_node *d, t_edge *b, phydbl *best_score, t_tree *tree)
792 {
793 int i;
794
795 if(d->tax) return;
796
797 for(i=0;i<3;i++)
798 {
799 if((d->v[i] != a) && (!d->v[i]->tax))
800 {
801 if((d->b[i]->nni->score > *best_score-1.E-10) &&
802 (d->b[i]->nni->score < *best_score+1.E-10)) /* ties */
803 {
804 d->b[i]->nni->score = *best_score+1.;
805 }
806
807 if(d->b[i]->nni->score < *best_score)
808 {
809 *best_score = d->b[i]->nni->score;
810 }
811 }
812 }
813 }
814
815 //////////////////////////////////////////////////////////////
816 //////////////////////////////////////////////////////////////
817 /*
818 v
819 |
820 |
821 a
822 / \
823 d u
824 / \
825 v1 v2
826 */
NNI_Traversal(t_node * a,t_node * d,t_node * v,t_edge * b,int opt_edges,t_tree * tree)827 void NNI_Traversal(t_node *a, t_node *d, t_node *v, t_edge *b, int opt_edges, t_tree *tree)
828 {
829 int i,dir1, dir2;
830
831 /* printf("\n. a: %d d: %d b->is_alive ? %d -- [%d - %d]", a->num,d->num,b->is_alive,a->tax,d->tax); */
832
833 if(d->tax == YES)
834 {
835 if(opt_edges == YES) Br_Len_Opt(&(b->l->v),b,tree);
836 return;
837 }
838 else if(a->tax == YES)
839 {
840 if(opt_edges == YES && a->tax == YES) Br_Len_Opt(&(b->l->v),b,tree);
841 for(i=0;i<3;++i)
842 if(d->v[i] != a)
843 {
844 Update_Partial_Lk(tree,d->b[i],d);
845 NNI_Traversal(d,d->v[i],a,d->b[i],opt_edges,tree);
846 }
847 Update_Partial_Lk(tree,b,d);
848 }
849 else
850 {
851 tree->n_edges_traversed++;
852 NNI_Core(a,d,v,b,opt_edges,tree);
853
854 dir1 = dir2 = -1;
855 for(i=0;i<3;++i)
856 if(d->v[i] != a)
857 {
858 if(dir1 < 0) dir1 = i;
859 else dir2 = i;
860 }
861
862 Update_Partial_Lk(tree,d->b[dir1],d);
863 NNI_Traversal(d,d->v[dir1],a,d->b[dir1],opt_edges,tree);
864 Update_Partial_Lk(tree,d->b[dir2],d);
865 NNI_Traversal(d,d->v[dir2],a,d->b[dir2],opt_edges,tree);
866
867 Update_Partial_Lk(tree,b,d);
868 }
869 }
870
871 //////////////////////////////////////////////////////////////
872 //////////////////////////////////////////////////////////////
873
874
NNI_Core(t_node * a,t_node * d,t_node * v,t_edge * b,int opt_edges,t_tree * tree)875 void NNI_Core(t_node *a, t_node *d, t_node *v, t_edge *b, int opt_edges, t_tree *tree)
876 {
877 phydbl lk0,lk1,lk2;
878 phydbl rnd;
879 t_node *v1,*v2,*u,*dum;
880 scalar_dbl *l0,*l1,*l2;
881 phydbl p_accept;
882 int i;
883
884
885 l0 = l1 = l2 = NULL;
886
887 lk0 = UNLIKELY;
888 lk1 = UNLIKELY;
889 lk2 = UNLIKELY;
890
891 v1 = v2 = NULL;
892 for(i=0;i<3;++i)
893 if(d->v[i] != a)
894 {
895 if(v1 == NULL) v1 = d->v[i];
896 else v2 = d->v[i];
897 }
898 assert(v1 != NULL);
899 assert(v2 != NULL);
900
901 dum = NULL;
902 rnd = Uni();
903 if(rnd < .5)
904 {
905 dum = v1;
906 v1 = v2;
907 v2 = dum;
908 }
909
910 u = NULL;
911 for(i=0;i<3;++i) if(a->v[i] != d && a->v[i] != v) { u = a->v[i]; break; }
912
913 if(opt_edges == YES) Br_Len_Opt(&(b->l->v),b,tree);
914 lk0 = Lk(b,tree);
915 l0 = Duplicate_Scalar_Dbl(b->l);
916
917 /* Swap_Partial_Lk_Extra(b,a,0,tree); */
918 /* Swap_Partial_Lk_Extra(b,d,1,tree); */
919
920 // First NNI
921 Swap(v1,d,a,u,tree);
922 // Update partial likelihood looking up
923 Update_Partial_Lk(tree,b,a);
924 // Update partial likelihood looking down
925 Update_Partial_Lk(tree,b,d);
926 // Evaluate likelihood
927 if(opt_edges == YES) Br_Len_Opt(&(b->l->v),b,tree);
928 lk1 = Lk(b,tree);
929 /* if(lk1 > lk0 + tree->mod->s_opt->min_diff_lk_move) */
930 /* { */
931 /* tree->fully_nni_opt = NO; */
932 /* return; */
933 /* } */
934 l1 = Duplicate_Scalar_Dbl(b->l);
935 // Unswap
936 Swap(u,d,a,v1,tree);
937
938
939
940 // Second NNI
941 Swap(v2,d,a,u,tree);
942 // Update partial likelihood looking up
943 Update_Partial_Lk(tree,b,a);
944 // Update partial likelihood looking down
945 Update_Partial_Lk(tree,b,d);
946 // Evaluate likelihood
947 if(opt_edges == YES) Br_Len_Opt(&(b->l->v),b,tree);
948 lk2 = Lk(b,tree);
949 /* if(lk2 > lk0 + tree->mod->s_opt->min_diff_lk_move) */
950 /* { */
951 /* tree->fully_nni_opt = NO; */
952 /* Free_Scalar_Dbl(l1); */
953 /* return; */
954 /* } */
955 l2 = Duplicate_Scalar_Dbl(b->l);
956 // Unswap
957 Swap(u,d,a,v2,tree);
958
959 /* Swap_Partial_Lk_Extra(b,a,0,tree); */
960 /* Swap_Partial_Lk_Extra(b,d,1,tree); */
961
962 /* if((u->tax == YES && !strcmp(u->name,"tax57")) || */
963 /* (v1->tax == YES && !strcmp(v1->name,"tax57")) || */
964 /* (v2->tax == YES && !strcmp(v2->name,"tax57"))) */
965 /* printf("\n. lk0: %G lk1: %G lk2: %G l0: %G l1: %G l2: %G",lk0,lk1,lk2,l0->v,l1->v,l2->v); */
966
967 /* printf("\n. a: %d d: %d -- lk0: %f lk1: %f lk2: %f p: %G %G %G %G",a->num,d->num,lk0,lk1,lk2,p_accept,l0->v,l1->v,l2->v); */
968
969 p_accept = exp((lk1-lk0)/(tree->annealing_temp+1.E-6));
970 if(Are_Equal(lk1,lk0,tree->mod->s_opt->min_diff_lk_local) && Are_Equal(tree->annealing_temp,0.0,1.E-3)) p_accept = .0;
971 rnd = Uni();
972
973
974 if(rnd < p_accept && lk2 < lk1)
975 {
976 Swap(v1,d,a,u,tree);
977 Copy_Scalar_Dbl(l1,b->l);
978 tree->c_lnL = lk1;
979 Update_Partial_Lk(tree,b,a);
980 Update_Partial_Lk(tree,b,d);
981 tree->fully_nni_opt = NO;
982 }
983 else
984 {
985 p_accept = exp((lk2-lk0)/(tree->annealing_temp+1.E-6));
986 if(Are_Equal(lk2,lk0,tree->mod->s_opt->min_diff_lk_local) && Are_Equal(tree->annealing_temp,0.0,1.E-3)) p_accept = .0;
987 rnd = Uni();
988 if(rnd < p_accept)
989 {
990 Swap(v2,d,a,u,tree);
991 Copy_Scalar_Dbl(l2,b->l);
992 tree->c_lnL = lk2;
993 Update_Partial_Lk(tree,b,a);
994 Update_Partial_Lk(tree,b,d);
995 tree->fully_nni_opt = NO;
996 }
997 else
998 {
999 Update_Partial_Lk(tree,b,a);
1000 Update_Partial_Lk(tree,b,d);
1001 Copy_Scalar_Dbl(l0,b->l);
1002 tree->c_lnL = lk0;
1003 }
1004 }
1005
1006 Update_PMat_At_Given_Edge(b,tree);
1007
1008 if(l0) Free_Scalar_Dbl(l0);
1009 if(l1) Free_Scalar_Dbl(l1);
1010 if(l2) Free_Scalar_Dbl(l2);
1011 }
1012