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